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ABSTRACT 



A quantitative analysis was carried out to determine the 
stresses present in the leading-edge segment of a P-3C 
aircraft operating within and outside the normal operating 
envelope of the aircraft. The purpose of the analysis was to 
ascertain whether a specific weakness may exist in the 
leading-edge structure which might endanger future operating 
flight crews. A three-step process consisting of a static 
aeroelastic span-load analysis, an inviscid two-dimensional 
panel method, and finite-element analysis was employed in the 
course of the evaluation. Lift-coefficient distributions from 
the wing span-load analyses were used in the two-dimensional 
panel method to determine the pressure distribution around the 
leading edge, which was then used as input to the finite- 
element analysis. Additionally, static aeroelastic-derived 
wing-twist effects were included in the structural model. The 
results of the analysis suggest that the leading edge segment 
studied may experience stress levels sufficient to cause 
failure within the normal operating envelope. 



iii 



ij 



THESIS DISCLAIMER 

Although certain commercial software products are 
referred to in the body of this thesis, such references do not 
constitute recommendations for use or endorsements of the 
products . 



iv 



DUDLEY KNOX LIBRARY 
rsLWAL POSTGRADUATE SCHOOL 
MONTEREY CA 93943-6101 



TABLE OF CONTENTS 



I. INTRODUCTION 1 

II. GENERAL AERODYNAMICS 4 

A. MISHAP AIRCRAFT SPECIFICATIONS 4 

B. PROBLEM STATEMENT 6 

C. ASSUMPTIONS 7 

III. WING SPAN-LOAD ANALYSIS 9 

A. THEORY 9 

B. USAGE OF THE SPAN-LOAD PROGRAM 15 

1. Additional Loading 16 

2. Built-in Geometric Twist 18 

3. Engine Dead-weight Twist 19 

4. Camber 20 

5. Aileron Float 21 

6. Aileron Deflection Angle 21 

, 7. Roll Helix Angle 21 

8. Propeller Slipstream Effect 22 

9. Tail Contribution 23 

C. EFFECTS NOT INCLUDED 25 

1. Propeller Gyroscopic Moment 25 

2. Wing Fuel Moment 26 

D. APPLICATION 26 



V 



IV. TWO-DIMENSIONAL PRESSURE ANALYSIS 30 

A. THEORY . 30 

1. Coordinate Axis and Panel Numbering System 31 

2. Flow Formulation 31 

3. Boundary Conditions 33 

4. Influence Coefficients 34 

5. Numerical Solution Method 35 

6. Velocity and Pressure Distribution .... 36 

B. APPLICATION 37 

1. P-3 Airfoil Section 37 

2. Program Inputs and Outputs 37 

3. Program Operation 41 

4. Verification of the Program 42 

5. Flight Regime Selection 42 

a. Method Employed 43 

b. Symmetric Airfoil 44 

c. Cambered P-3 Airfoil 44 

d. Analysis 45 

V. FINITE ELEMENT ANALYSIS 46 

A. THEORY 46 

1. Accuracy 47 

2. Equations 48 

B. STRUCTURAL REPRESENTATION 49 

1. The Leading Edge Structure 49 

2. The Model 51 



VI 



C. APPLICATION 



57 



1. Input Data 58 

2. Finite Element Analysis Results 59 

a. Load Case 1 59 

b. Load Case 2 62 

c. Load Case 3 71 

d. Load Case 4 72 

e. Load Case 5 73 

f. Load Case 6 73 

3. Discussion 74 

a. Extending the Structure 76 

b. Stress Concentrations 76 

VI. CONCLUSIONS AND RECOMMENDATIONS 79 

A. CONCLUSIONS 79 

B. RECOMMENDATIONS 81 

1. Validation 81 

2. Flight Envelope Modifications 82 

LIST OF REFERENCES 84 

APPENDIX A SPAN-LOAD ANALYSIS PROGRAM 86 

APPENDIX B AIRFOIL PANEL METHOD PROGRAM 98 

INITIAL DISTRIBUTION LIST 108 



vii 



ACKNOWLEDGEMENT 



This thesis would not have been possible without the help 
and support of many. Foremost among these are the members of 
my family, Fran, Casey and Ryan, who have endured my absences, 
aggravation and absent-mindedness. Your love has sustained 
me. 

To the scholar from who's seemingly endless reserve of 
talent and patience I have drawn for many months, I can only 
say, "Thank you. Professor Lou Schmidt". 

Thanks also to those who have provided help, inspiration, 
equipment and time in this endeavor; LCDR Rick Bobbitt who 
brought this topic to my attention and gave of his time and 
advice; Dr. Hank Smith who co-sponsored Rick's proposal; 
Professor Rick Howard who helped with the simulator run and 
photography, as well as scholarly advice; Professor Ed Wu for 
his inspiration and aid in the area of finite element 
analysis; Professor Gerald Lindsey for his help in structural 
analysis; and Nam Phan and Li Chen, who provided the 
foundation and support for the finite element model. 



viii 



I. 



INTRODUCTION 



On 13 February 1988, a U.S. Navy (USN) P-3C Update I Orion 
aircraft experienced an in-flight failure of a wing leading- 
edge section during a high-speed, low-altitude maneuver. The 
aircraft was able (with some difficulty) to return to its 
departure point and make a safe landing. Some three years 
later, on 27 April 1991, an Orion operated by the Royal 
Australian Air Force (RAAF) suffered a similar but more 
extensive in-flight failure when it lost three wing leading- 
edge sections during a maneuver similar to the U.S. Navy 
mishap. The Australian crew was unable to return to the 
runway due to the loss of lift available from the damaged 
wings. The aircraft impacted the ocean surface short of its 
island destination in a nose high-attitude with maximum power 
set on all four engines. One crewmember was killed when a 
propeller tore loose from its engine and entered the fuselage 
where he was sitting. 

The leading edge is the rounded front portion of the wing 
which is physically attached to the front wing spar and is 
essential in the production of lift by an airfoil. Loss of a 
.leading-edge segment results in a dramatic reduction in wing 
lift capability combined with a corresponding increase in 



1 



drag. The P-3 has three of these leading-edge segments on 
each wing, separated by the engine nacelles. They are best 
referred to as inboard, center and outboard sections. Figure 
1 shows the P-3 aircraft. The inboard sections are those 
between the fuselage and the 




sections cover the remaining distance from the outboard 
engines to the wing tips. The length (fore and aft) of these 
leading-edge sections is 15% of the chord (total fore and aft 
distance) of the wings. 

The USN P-3 lost the starboard wing center leading-edge 
section while the RAAF aircraft lost both of its center 
sections and the starboard inboard section. This study 
focuses on the center leading edge sections. 

While both of these mishaps were investigated by the 
appropriate authorities, the cause of the failures were 
undetermined, though widely suspected to be the result of 
aircraft overstress due to pilot error. In both cases, the 
pilots did not feel that they had exceeded the normal 
operating envelope (this envelope will be detailed later) for 
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the aircraft. Even if the aircraft were operated outside the 
authorized envelope, the question remains as to why the 
leading edge sections failed before some other component. 

This thesis was undertaken for the purpose of studying the 
aerodynamic loading and structural response of a leading edge 
section due to operation within and outside the normal flight 
envelope. Ultimately, its aim was to determine whether there 
may be a need to further restrict the operating envelope or 
recommend some modification to the existing structure in order 
to prevent a further recurrence of the in-flight failure. 
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II. GENERAL AERODYNAMICS 



A. MISHAP AIRCRAFT SPECIFICATIONS 

The configurations and operational parameters of the 
mishap aircraft were similar in that both were at high gross 
weight (RAAF at 127,000 pounds, USN at 135,000 pounds) and 
operating at high airspeeds (RAAF approximately 380 knots, USN 
approximately 350 knots). Both aircraft executed a pullup 
maneuver from an altitude of less than 400 feet at the time of 
leading edge failure. The RAAF maneuver was a straight pullup 
(wings level), while the USN maneuver was a starboard rolling 
pullup (meaning that the aircraft was banked to the right as 
the pullup maneuver was executed) . Interviews with some USN 
crewmembers indicate that the failure of their leading-edge 
section may have begun a few seconds earlier as the aircraft 
rolled from a right bank to wings level while inbound for the 
above stated pullup. This, they said, was reported to them by 
ground observers. Table 1 presents the relevant P-3C 
dimensional data used during this analysis while Table 2 
delineates the aircraft performance parameters. Figure 2 
shows the operating flight envelope for the aircraft in a 
flaps-up, gear-up configuration. 
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TABLE 1. 



DIMENSIONAL DATA FRefs. 1,2,31 



Wing 

Area, S (ft^) 1300 

Span, b (theoretical, ft) 99 

MAC, cbar (ft) 14.1 

Aspect Ratio, A 7.5 

Taper Ratio, X 0.4 

Dihedral, (.25c„, degrees) 5.0 

Incidence, Root (degrees) 3.0 

Tip 0 . 5 

Airfoil Section, Root NACA0014-1 . 10 40/1.051 cl^=.3,a=.8 

Tip NACA0014-1. 10 40/1.051 cli=.4,a=.8 
Straight Element 0.15c 

Chord, Root (ft) 18.9 

Tip 7 . 6 

Aileron 

Area, S^ (ft^) 45.5 

Hinge Line (c„) 0.725 

Deflection Limit, Up (degrees) -23.3 

Down +16.2 

Horizontal Tail 

Area, ( ft^ ) 321.8 

Tail Length (.25cbar„ to .25cbar^, ft) 49.8 



TABLE 2. PERFORMANCE PARAMETERS FRefs. 1.2.31 



Flight Design Gross Weight (pounds) 

Load Factor (through 135,000 pounds, G) 
Maximum Operating Speed, Sea Level (knots) 
Center of Gravity Limits (135,000 lbs, %MAC) 
Maximum Shaft Horsepower (per engine) 

Maximum Lift Coefficient, Cn^ax (po^/er off) 
Lift Curve Slope, C^^/ Tail Off (power off) 

Tail On (power off) 

Moment Slope, C^^l/ Tail Off, c.g @ 0.25cbar^^ 



135000 

-1.0 to +3.0 
405 

21.5 to 31.0 
4600 
1.30 
4.84 
5.50 
0.19 
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[Ref c 2 ] , 



B. PROBLEM STATEMENT 

Accurate analysis of the problem within limited time and 
budget constraints available for completion of the work served 
to restrict the options for solution methods. While 
instrumented flight and laboratory tests would probably be 
more precise, these choices were not feasible. It was 
necessary to develop an analytic method which could provide 
credible results within the academic environment. 

Determination of the loads applied and their effect upon 
the leading edge structure was accomplished by a three-step 
process. The first step was to conduct a static aeroelastic 
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span-load analysis of the wing using a computer program to 
determine the section lift coefficients and structural twist 
on the P-3 wing. Next, a two-dimensional panel method was 
employed to find the pressure distribution around the leading 
edge. Finally, the forces derived from the pressure 
distribution were used as the loads applied in a finite 
element analysis computer application. In addition, 
consideration was given to the aeroelastically-derived 
spanwise wing twists which were introduced as torsion-like 
deflections in the finite element analysis. 

C. ASSUMPTIONS 

Certain assumptions were made in conducting the analysis. 
While static aeroelasticity effects upon the wing were 
accounted for, the fuselage and tail were assumed as rigid 
structures. In addition, the effect of fuselage interference 
on wing lift distribution was neglected. That is, the wing 
was taken to be fully effective in producing lift even in the 
central region where it is influenced by the fuselage. The 
unswept, straight-tapered wing with an aspect ratio of 7.5 was 
modeled structurally in the span-load analysis by an elastic 
axis. Chordwise bending of the wing was neglected. Inviscid 
solution methods were employed in the span-load and airfoil 
analysis programs. It was felt that these assumptions would 
model the flowfield without inducing an unacceptable level of 
error in the overall outcome, and that this analysis was a 
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preferable method of solution to performing a computational - 
fluid-dynamics analysis where static aeroelastic influence 
could not be included. Static aeroelastic effects upon wing 
loadings were suspected of playing a large role in the problem 
because of the location of the wing's elastic axis at a 
constant 40 percent of chord, according to available 
information. This is well aft of the 25 percent of chord 
location at which a majority of the aerodynamic loadings are 
presumed to act, creating the potential for a significant 
coupling influence by structural deflection during the 
development of lift. 
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III. WING SPAN-LOAD ANALYSIS 



A. THEORY 

Because the P-3 wing (as well as nearly all others) is 
flexible, the aerodynamic- load distribution varies from that 
which would be seen when considering the wing as a rigid 
structure. Allowing for structural deflection of the wing 
increases the accuracy of the result in the attempt to model 
the true physics of the problem. The computer program used 
for this analysis, written in Microsoft Quick BASIC®, was 
based upon the work of Schmidt [Ref. 4]; c.f. Appendix A. The 
basic concepts employed in the program are presented below. 

Defining a set of linear, simultaneous equations for the 
span-load solution on a wing at a specified angle of attack 
starts with the following: 

+ ai2(^/q)2 +••• + ain(^/q)n = 

where the wing is cut into a series of spanwise stations and 

• a^^j = aerodynamic influence coefficient of (i/q) at 

station "j" upon induced angle at control station 

II ^11 

• £ = running span-load (lb-in“^) 

• q = freestream dynamic pressure (Ib-in"^) 

• ai = geometric angle of attack at control station "i". 
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The relationship states that the span-load-induced downwash 
velocities satisfy flow tangency at a control point. 
Expressed as a matrix equation this system becomes: 

[A]{i/q} = {a} 



where , 



• [A] = square (n x n) matrix of aerodynamic influence 

coefficients ^ength"^) 

• = column (n x 1) matrix of span-load values 

(length) 

• {a} = column (n x 1) matrix of angle-of-attack input 

(radian) 



In similar fashion, a matrix equation relating the 
structural twist at a wing station "i" to the span-loading may 
be developed from: 



+. . . + 

or, q(Sij(je/q )2 + s^ 2 (je/q)j +... + s^„(je/q)„) = 

which may be stated to apply to all the wing stations as 

before to yield the matrix equation, 



where , 



q[S]{i/q} = { 2 ^ag} 



• [S] = square (n x n) matrix of structural influence 

coefficients (length-force"^) 

• {AXg} = column (n x 1) matrix of effective angle-of-attack 

changes due to structural twist (radian) 



Next, the two matrix equations may be combined to generate 
a single equation for an elastic wing as follows: 
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(A){-e/q)E = {“>R +{^)s 
[A)(^/q)E = {a)R + q[S]{i/q)E 
[ [A] - q[S]l {je/q)E = {a}R 

yielding the span-load solution for the elastic wing, 

{^/q)E = [ [A] - q[S]] -l{a)R 

where the subscript "E" is used to mean elastic and "R" to 
mean rigid. 

The concept of mathematical symmetry may be employed to 
develop equations for a symmetric and anti- symmetric load case 
by considering the wing to be divided into left and right hand 
wing panels at the spanwise centerline. This process allows 
a superposition of linear aerodynamic solutions such as the 
combination of a symmetric pullup and anti-symmetric roll 
rate to yield the total solution for a rolling pullup. It 
also reduces computing time since the size of the matrices is 
half the original required for the total wing. In employing 
this technique it is necessary to distinguish between the 
symmetric and anti-symmetric forms of the aerodynamic 
influence coefficients. 

The approach employed for developing the aerodynamic 
influence coefficients involves use of a technique known as 
the "Modified Weissinger" approach, wherein the wing is 
divided into a series of spanwise stations with swept bound 
vortices attached at the local quarter-chord point, giving 
rise to horseshoe vortices extending downstream to infinity in 
accordance with Helmholtz' laws. The vortex strengths are 
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determined in accordance with the Biot-Savart law, such that 
the induced downwash angle at all local three-quarter chord 
points (from the summation of all vortices) are equivalent to 
the local geometric angle of attack. The three-quarter chord 
points are termed "control points". Enforcement of the stated 
boundary conditions replicates the occurrence of tangential 
flow over the surface of the wing. Application of the Biot- 
Savart and Kutta-Jukowski laws to the geometric relationships 
of the horseshoe vortices and control points results in 
development of the symmetric and anti -symmetric aerodynamic 
influence coefficient matrices, 

[ -^s ] LH^ 

where the subscripts "RH" and "LH" denote the right and left 
hand wing panels, respectively. These influence coefficient 
matrices may next be substituted in the previously-developed 
equations to yield, 

[Ag]{i/q} = {Og} 

[Ag]-{i/q} = {ag} 

Subsonic compressibility effects were included in the 
development of the aerodynamic influence coefficients using 
the Prandtl-Glauert planform distortion approach [Ref. 5]. In 
this method, the chordwise dimension of the planform is 
increased according to the relationship. 
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X 






v/l-W‘ 



The "Modified Weissenger" approach was altered to a panel 
form for this analysis in that the wing half was divided into 
five chordwise and ten spanwise stations, requiring 
manipulation of (50 x 
50) matrices. A 

representative sketch of 
the method employed is 
shown in Figure 3. The 
vector U in the figure 
represents the free 
stream velocity while r 
is the circulation 
strength of the vortex 
filament . 

Development of the 
structural influence 
coefficient matrix in 
the program is accomplished through determination of the 
moments exerted about the longitudinal (X) and lateral (Y) 
axes of the wing at points along the wing elastic axis as 
shown in Figure 4. These points correspond to the mid-span 
locations of the individual panels, along the elastic axis. 
These moments are found by solving the matrix equations: 
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(My)=q(B/N)[t]{je/q) 

(M^)=q(B/N)2[mHje/q) 

Where B is the span, N 
is the number of 
spanwise local stations, 

[t] is a torsional 
matrix and [m] is a 
bending moment matrix. 

From these equations, a 
coordinate axis rotation can be performed as follows: 

{T} = cosi^{My} + sini^{M^} 

{M} = - sini^{My} + cosi^{Mjj} 

Next, torsional and bending stiffness (GJ and El) values along 
the elastic axis are employed to determine the angular 
deflections resulting from the torsion and bending moments 
from the equations. 




y 

COSAg 



0t(y)= 



I 



T(z) 

GJ{z) 



dz 



y 

cosX. 



©t(/n) = 



l 



M{z) 

Eliz) 



dz 



These are then discretized as. 






B 

2WcosAg 



N 

( 

J=I + 1 




M{J) ^ M{I) 
EI[J) EI{I) 



A0t(I) 



B 

2NcosAg 



N 

( E 

J=J+1 



T(J) ^ T(I) 
GJ(J) GJ(I) 
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or expressed in matrix form, 






2 iVcosA^ ^ 



The angle-of-attack column vector may then be found as, 

{Axg} = cosi^{0t} - sini^{0jjj} 

which is also: {^ 3 } = q[S]{l/q}. 



The preceding equations may then be tied together to form the 
equation for the structural influence matrix, 

[S] = (B/2N)2[uJ [1/GJ][ cosi^[t] + sin^(B/2N) [m]] 

+ tani^(B/2N)2[u] [1/EI][ sini^[t] - cosi^(B/2N) [m]] 

B. USAGE OF THE SPAN-LOAD PROGRAM 

Application of the wing span-load program required the 
determination of the following influences: 

• Additional loading distribution due to wing angle of 
attack 

• Built-in geometric twist 

• Airfoil camber distribution 

• Dead-weight induced wing twists due to propulsion system 
weight 

• Aileron float angle 

• Propeller slipstream effects 

• Aileron control deflection 

• Roll helix angle 
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The last two influences involved anti -symmetric wing span-load 
solutions. The total wing span-load distribution, which 
provided a measure of wing lift coefficient (C^)/ moment 
coefficient and section lift coefficient (Cjg), was 
obtained by an appropriate linear combination of the above 
influences. These influences were incorporated in the program 
as discrete angle of attack adjustments at the control points. 
Finally, a specified airplane lift coefficient required an 
estimate of the tail lift contribution before the 
representative wing Cl could be estimated. The tail lift 
contribution was based upon trimming the P-3 airplane for a 
specified flight condition and the assumption that the body 
and horizontal tail behaved approximately as rigid structures. 

1 . Additional Loading 

Static aeroelastic effects upon wing span-loading due 
to geometric angle of attack were incorporated in the program 
as a selectable input from the operator at the computer 
terminal. Calculation of the appropriate angle of attack for 
a given flight condition was based on the fundamental equation 

with appropriate modification for tail lift contribution as 
delineated in subsection nine of this Chapter. 

Early in the analysis, it was found that static 
aeroelastic effects had a dramatic impact on additional 
loading. Solutions were obtained using the span-load program 
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for a wide range of dynamic pressures. The q equals zero 
solution corresponded to a rigid-wing case. A dynamic 
pressure of 3.7 psi corresponded too the aircraft operating at 
a Mach number (M) of 0.6 at sea level. The variation of wing- 
alone Clq with dynamic pressure, shown in Figure 5, indicates 
a 41 percent increase due to static aeroelastic influences at 
the sea level flight condition. The increase in lift-curve 
slope is associated with a spanwise variation of structural 
twist as shown in Figure 6. At a Mach number of 0.6 for sea- 
level flight, each degree of geometric angle-of-attack input 
at the wing root results in 1.95 degrees of a at the tip, with 
the added 0.95 degrees being due to the structural twist 
component in the a direction. The static aeroelastic 
influence upon the wing aerodynamic center was determined as 



T3 

O 

w. 

w. 

(V 

Q. 



O 

(V 

Q. 

JO 

in 

(V 

3 

O 




Dynamic Pressure, 0 (psi) 



Figure 5. Aeroelastic variation of lift-curve slope 
[Ref. 6] 
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being negligible, a result which may be attributed to the wing 
0.25 chord line having a small sweep angle. 




Fraction of Wing Semispan, 2y/B 

Figure 6. Aeroelastic variation of twist 
[Ref. 6] 

2. Built-in Geometric Twist 

The effect of 2.5 degrees of washout was included in 
the program by linearly varying the local (panel) angle of 
attack moving outward from zero at the root to -2.5 degrees at 
the tip. Figure 7 depicts the effect of washout on section 
lift and twist distribution for the rigid and elastic P-3 
wing cases. The magnitude of the negative section Cjg values 
is seen, in Figure 7, to increase by 21 percent while the wing 
tip experiences a 1.5-degree negative twist due to aeroelastic 
influences . 
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Figure 7. Effect of built-in geometric twist on section 
lift coefficient and twist distribution. 



3. Engine Dead-weight Twist 

A separate version of the span-load program was 
developed and run to determine the wing twist induced per G 
due to the effect of engine dead-weight moment. Output 
information was generated in terms of discretized angle-of- 
attack adjustments and included in the basic span-load 
program. An engine and propeller assembly weight of 3974 
pounds was assumed to act at the (x,y) coordinates (coordinate 
system as shown in Figures 3 and 4) of (-51", 187") and (-45", 
357") for the inboard and outboard engines, respectively. 
Figure 8 shows the effect of engine dead-weight twist on the 
elastic twist distribution at three G's and 405 knots. The 
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wing is twisted to a maximum of -4.5 degrees at the tip under 
this load condition. 




Figure 8. Effect of engine dead-weight on structural 
twist . 

4 . Camber 

Although listed in Table 1 as a NACA 0014 which tapers 
to a NACA 0012, the airfoil used for the P-3 wing is not a 
symmetric section, as indicated by available data [Ref. 8, p. 
2-470]. It is, in fact, a hybrid with considerable camber. 
The wing camber has an effect on aeroelastic behavior and was 
included in the analysis by determination of camber-line slope 
at the five chordwise control points which make up the wing 
model and entering the negative of this value as an adjustment 
to the station angle of attack. 
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5. Aileron Float 



The span-load analysis included the effect of aileron 
float as given in a graph of aileron angle versus airspeed in 
a Lockheed report [Ref. 9 , p. 123]. An equation curve fit for 
the 50-pound control-wheel-force curve was developed and 
incorporated in the program. Deflected aileron surface area 
was matched by a selection of panels on the wing model for 
varying degrees of deflection in order to closely emulate the 
aircraft control-surface deflection and aeroelastic effect. 

6. Aileron Deflection Angle 

Using the same process as that used for the aileron 
float angle, the effect of aileron deflection for 
consideration of anti-symmetric solutions to emulate roll 
maneuvering was also included. A curve fit to the data given 
for available 50-pound control-wheel-force deflection in the 
same report [Ref. 9 , p. 123] was used to generate the 
deflection angles. 

7. Roll Helix Angle 

Data [Ref. 9 , p. 118] for available tip helix-angle 
(pb/2V) variation with airspeed were curve fitted and 
incorporated as another control point angle-of-attack 
variation in the program. A roll helix angle of 2.63 degrees 
provided roll moment equilibrium with the available 9.3 
degrees of aileron deflection at 275 knots equivalent 
airspeed, at a control wheel force of 50 pounds. At 405 



21 



knots, the roll helix angle was 1.03 degrees for roll 
equilibrium with an aileron deflection of 8.2 degrees. The 
merging of the anti-symmetric input due to aileron deflection 
with the symmetric pullup solution allowed a comprehensive 
analysis of the static aeroelastic effect of a rolling pullout 
maneuver as described in part A of this Chapter. An example 
of the individual and combined effects of these components on 
wing section lift coefficient { 0 ^) distribution is shown in 
Figure 9, where the outer corner of the operating envelope for 
roll maneuvering (405 knots, 2.4 G's) was examined. 



Rolling Pullup Contributions 




Wing Station (2y/b) 



Figure 9. Rolling pullup contributions to section lift 
distribution . 



8. Propeller Slipstream Effect 

The increase in dynamic pressure over the wing due to 
the propellers was calculated using momentum theory as stated 
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in Glauert [Ref. 10, p. 200] according to the following 
formula: 

T = Ap(V + Vi)V2 

where , 

• T = propeller thrust (pounds) 

• A = propeller disk area (ft^) 

• V = freestream velocity (ft/sec) 

• = velocity increase behind the propeller, determined 

from known thrust. 

The velocity increase, v^^, was converted to a dynamic pressure 
boost and incorporated in the span-load program in the area 
behind the propellers. 

9. Tail Contribution 

Horizontal tail and fuselage moment effects on wing 
span-load distribution resulted in an adjustment to the wing 
angle-of-attack value used as input to the program. This 
adjustment was calculated according to the relationship, 

“ ^Ladd ^Lt ^LO 

where, 

• Cl = lift coefficient required for flight condition 

• CLadd “ lift coefficient due to additional loading (due to 

angle of attack) 

• CLt = tail contribution to lift coefficient 

• Clo = tail-off lift coefficient at zero angle of attack 

(due to camber, wing twist and dead-weight twist) 
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CLt = -AC^t(cbar/£t) 



and, 
where , 

• = tail-moment coefficient 

• cbar = mean aerodynamic chord of the wing 

• = distance from .25cbar„ to .25cbar^ 

Additionally, at airframe trim = 0), 

^mt ~ *^m0 *“mCL^Ladd 

where, 

• Cjj^Q = moment coefficient at zero angle of attack. 

• = tail-off airplane dCj^^/dC^, c.g. = 0.25cbar 



The above was assembled and solved for C^add give. 



-Ladd' 



■‘■t 

Z~~charZ 

■‘•^—2 — ^mCL 



The value of tail-off bhe aircraft was found from 

wind tunnel data to be approximately 0.19 [Ref. 3, p. 20]. 
The value of C^^dd given flight condition was then used 

to solve for the angle of attack in the span-load program 
according to the formula: 

® “ ^Ladd/^La 
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The value of was available as an output, for the elastic 
wing, from the program for any given airspeed by entering one 
radian for the angle of attack input. 

A correction of +0.03 to C^j^q due to fuselage effect 
was made after finding that the value in the wind tunnel 
data and that from the program differed by this amount. This 
correction was verified by calculations made in accordance 
with Etkin [Ref. 7, p. 334] concerning the effect of body and 
engine nacelles on neutral point location. 

C. EFFECTS NOT INCLUDED 

Two other factors were considered as possible contributors 
to the static aeroelastic problem, but were found to have no 
significant impact and therefore not included in the span-load 
program. These factors were propeller gyroscopic effects and 
the effect of moments due to wing fuel. 

1. Propeller Gyroscopic Moment 

This phenomenon was investigated using the following 
formulation from [Ref. 11]: 

M = 2IpOJfi 
Ip = mk^ 

where, 

• M = moment due to gyroscopic precession (Ib-ft) 

• Ip = polar moment of inertia of each blade (Ib-ft sec^) 

• CO = rotation rate of the propeller (rad-sec“^) 
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• fi = pitch rate of the aircraft (rad-sec“^) 

• m = mass of propeller blade (slugs) 

• k = radius of gyration (ft) 

The analysis showed that, allowing for a f ive-degree-per- 
second aircraft pitch rate, the moment developed by each 
propeller was 1225 Ib-ft in a counterclockwise direction as 
viewed from above the propeller. This moment was considered 
insignificant since it is applied in the lateral plane and 
does not influence the static aeroelastic span-load solution. 

2e King Fuel Moment 

Analysis of the fuel tank geometry and location 
revealed that the center of mass of the fuel (with full wing 
tanks) lies nearly coincident with the elastic axis. Any 
torsional moment derived from this source would be negligible. 

D. APPLICATION 

Linear summation of the contributing factors addressed 
earlier in this Chapter resulted in a prediction of the total 
span-load distribution for the P-3 wing under any given flight 
condition. The primary flight conditions of concern in this 
analysis were those encountered during the aircraft mishaps as 
discussed in the introduction. The basic premise applied was 
to examine limit loads at the edge of the operating envelope 
and then expand the analysis to regions outside the envelope, 
at the ultimate load condition. In addition it was decided to 
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first examine the loads encountered in a symmetric pullup of 
three G's at 275 knots. This speed was chosen because of an 
observed higher load on the leading edge than determined at 
the upper right corner of the flight envelope, at 405 knots. 
(A further discussion of this loading phenomenon will be 
presented in Chapter IV. ) Following that, the target flight 
condition was extended to 4.5 G's at 325 knots, which 
represented an approximate extension of the envelope using a 
value for 1.3. Next, the effect of anti -symmetric 
loading in the form of a starboard and then a port rolling 
pullup were considered in order to assess the contributions of 
aileron deflection and roll helix angle. Rolling pullup 
analyses were done at 2.4 G's and 275 knots to remain inside 
the envelope and maintain some congruity with the symmetric 
loading case. Presented in this section are plots of the 
impact of some of the various flight maneuvers on section lift 
coefficient and twist distribution, using the fully-developed, 
tail-on solution with all contributing factors included. 
Figures 10 and 11 show the 275-knot, 3-G symmetric condition, 
and compare section lift coefficient and structural twist for 
the rigid-wing and elastic-wing cases. In Figures 12 and 13, 
the same distributions are depicted for the 275-knot, 2.4-G 
rolling pullup load condition. 



27 




Figure 10. Lift coefficient distribution at 275 knots, 
3-G, symmetric pullup. 




Figure 11. Twist distribution at 275 knots, 3-G, 
symmetric pullup. 
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Figure 12. Lift coefficient distribution at 275 knots, 
2.4-G, rolling pullup. 




Figure 13. Twist distribution at 275 knots, 2.4-G, 
rolling pullup. 
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IV. TWO-DIMENSIONAL PRESSURE ANALYSIS 



A. THEORY 

Like the static aeroelastic span-load analysis, the two- 
dimensional pressure analysis employed in this thesis involves 
the application of linear superposition. This is a 
consequence of the pressure analysis being based upon a 
solution to the LaPlace equation, a linear, homogeneous 
second-order partial differential equation. Linearity allows 
the problem to be subdivided into three separate elements 
(which will be described later in this Chapter) and added 
together. The actual application was based upon a panel 
method similar to the technique used in the preceding Chapter 
except that now, instead of dividing a wing planform into 
chordwise panels, a two dimensional (x,y) airfoil is divided 
into panels along the perimeter of its surface. Camber is 
defined in the airfoil shape, instead of being added on as a 
discretized variation of angle of attack, as before. The 
analysis now concerns a vertical plane or cross-section. 
Steady (no variation in the flow field with time), inviscid, 
incompressible flow is assumed to exist. The panel method and 
formulation employed is documented in a Naval Postgraduate 
School thesis by Teng [Ref. 12]. 
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1. Coordinate Axis and Panel Numbering System 



The airfoil is considered to be fixed in an (x,y) 
coordinate system with its origin at the point of intersection 
of the chord line with the leading edge. The positive x axis 
points aft to the trailing edge while positive y is up. The 
panels which make up the airfoil surface are of varying 
lengths, depending mostly on the radius of curvature, and are 
numbered from 1 through n starting at the lower surface of the 
trailing edge and proceeding clockwise to the upper surface at 
the trailing edge. Delineating the end points of these panels 
are nodes which begin with the number 1 at the trailing edge 
and proceed along the same numbering path as the panels. The 
trailing edge point is counted twice, giving n+1 nodes in all. 

2 . Flow Formulation 

Consider some panel j on the surface of the airfoil . 
On this surface there exists a pair of singularity 
distributions, known as a source distribution qj and a 
vorticity distribution y. The strength of the source 
distribution varies from panel to panel while the vortex 
strength is the same for all panels. These singularity 
distributions satisfy LaPlace's equation and the far field 
boundary condition. 

Applying superposition, the overall flow field is 
considered to be made up of three individual flows and is 
represented by the equation. 
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^ = 0CO + 0s + <^>V 



where <p^ is the potential of the freestream flow, 

<f)^ = V„(x cosa + y sina) 

(pQ is the velocity potential of the source distribution of 
strength q(s) per unit length (s) and is calculated by. 



where (r) is the radial distance from some point at which a 
source and vortex flow exist, to the midpoint of the panel in 
consideration. In addition, <p^ is the velocity potential of 
a vorticity distribution of strength A(s) per unit length and 
is given by. 



where 0 is the angle formed by a line drawn along the radial 
distance (r) and the panel in question. 



the straight line which makes up each of the panels, where qj 
and Xare constant. The individual effects are then summed to 
give the total effect of the sources and vortices from all 
panels, as given in the equation. 





Each of the preceding equations is integrated along 




J~^ panel ( j) 



2n 2n 
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The calculation of $ requires solution of the (n+1) 
unknowns, qj (j = 1,2, ...,n) and y* This is accomplished 

numerically in the computer program. Once $ is known, the 

velocity can be found by taking the gradient (V) of The 

total velocity vector is found as, 

Vtotal = + V (</>3 + 0^) 

Next, the coefficient of pressure is found from the 
Bernoulli equation in the incompressible form. 



c„=i-( 



total 



2 



3 . Boundary Conditions 

Both the condition of flow tangency at the surface and 
the Kutta trailing edge condition must be satisfied as 
boundary conditions. As in the span-load analysis, control 
points are designated at which flow tangency is satisfied, 
except here the control point is taken as the mid point of the 
panel. It is stipulated that each control point will have 
tangential velocity, but that all normal velocities, 

will be exactly zero. 

The Kutta condition requires that the pressures on the 
upper and lower surface at the trailing edge be equal. Using 
Bernoulli's equation for steady potential flow, this state of 
pressure equilibrium is found to exist when the tangential 
velocities in the downstream direction are equal at the upper 



33 



and lower trailing-edge panels. In equation form this is 
written, 

The task then becomes one of using the boundary conditions to 
solve for the (n+1) unknowns. 

4. Influence Coefficients 

The concept of influence coefficients is again 
employed in this portion of the analysis, as it was in the 
span-load analysis. Here, the influence coefficients take 
the form of induced normal and tangential velocities at the 
control point of a given panel. These velocities are induced 
by the source and vorticity distributions of the other panels, 
and it is from this influence that they receive their 
designations: 

• = normal velocity induced at the i^^ panel control 
point by the source distribution on the panel. 

• A^^j = tangential velocity induced at the i^^ panel control 
point by the source distribution on the panel. 

• = normal velocity induced at the i^^ panel control 
point by the vorticity distribution on the panel. 

• = tangential velocity induced at the i^^ panel control 
point by the vorticity distribution on the panel. 

These influence coefficients are calculated through 

application of the geometric relationships which exist between 
the panels in conjunction with the formulas for the source and 
vorticity velocity potentials as given above. 
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5. Numerical Solution Method 



Using the influence coefficients, the boundary 
conditions can now be employed to write (n+1) equations which 
may be solved in matrix form for the (n+1) unknowns. The set 
of n equations comes from enforcement of the flow tangency 
boundary condition in the form, 

n n 

£ B"i^+V,sin(a-6i)=0 

Next may be written the enforcement of the Kutta boundary 
condition as, 

n n n n 

J=1 J=1 J=1 J=1 

The negative signs on the left side of the equation are due to 
the defined orientation of the tangential velocities as 
positive in the downstream direction. 

These equations may then be expressed in matrix form 
with the (n+1) unknowns (i.e., qj ( j=l , 2 , . . . , n) and y) 
arranged as a column (nxl) matrix multiplied by the 
(n+1 X n+1) influence coefficient matrix and set equal to a 
Bn+i column matrix. From this point, a Gaussian Elimination 
numerical technique may be employed to solve for the (n+1) 
unknowns [Ref. 13]. 
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6. Velocity and Pressure Distribution 

Once the qj and y found, the tangential velocities 

at the control points can be solved for according to the 
equation: 

n n 

Vtotai’E B%+v.cos(a-ei) 

J=1 J=1 

where i = 1,2, ...,n. From this, the individual pressure 
coefficients may be solved using the equation, 

(Cp)i = 1 - i=l,2,...,n 

At this point, the forces at work on the airfoil may be found 
by first integrating forces in the airfoil coordinate system 
as follows, 

n 

i=l 

n 

i=l 

Performing a coordinate axis rotation to re-align with that of 
the freestream yields the lift coefficient, 

Cjg = Cy cosa ” Cjj sina 

For this application, the values of 0,^ an Cy are found in 
discretized form at each panel as i ^ “ 

l,2,...,n, and then multiplied by the chord length and 
freestream dynamic pressure to arrive at the normal and 
chordwise force exerted at each panel. From this, the 
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leading-edge panels are selected and their forces collected 
for application as point loads in the finite element analysis 
portion. 

B. APPLICATION 

1. P—3 Airfoil Section 

As mentioned in Chapter 111 , the airfoil shape used in 
this application is as delineated by [Ref. 8]. Surface 
coordinate locations were solved using the tables and 
equations provided therein relative to a "wing reference 
plane" which appeared to correspond to a water line (i.e. the 
angle of incidence at the root was included in the 
definition). These coordinates were then rotated to an (x,y) 
coordinate system aligned with the chord of the airfoil as 
required in this panel method. The result was an airfoil 
consisting of some 48 panels, which was increased to 52 panels 
(53 node points) due to observed roughness of the leading edge 
shape when plotted. This smoothing of the leading edge shape 
was achieved by applying the specified leading edge radius to 
create intermediate node points. The basic shape of the 
airfoil, though not precisely to scale, is depicted in Figure 
14. The locations of the node points are also shown. 

2. Program Inputs and Outputs 

The desired output from the two-dimensional panel 
method program was a collection of forces and boundary 
constraints to be applied at node locations in the finite 
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Figure 14. P-3 airfoil showing panel nodes. (Not to 
scale ) 



element analysis phase. In order to achieve this, the program 
was altered considerably from the initial state as outlined 
earlier in this Chapter. The final FORTRAN code is available 
in Appendix B. Since it was necessary to develop loads to be 
applied to a three-dimensional model, the program was set up 
to iteratively compute the force distribution at a series of 
two-dimensional airfoil sections which ranged in location from 
the inboard end of the finite element model at wing station 
(WS) 256 to the outboard end at WS 320. (These locations 
correspond to the outboard 64 inches of the wing center 
section leading edge, located between the nacelles.) This 
section force distribution was scaled up from a chord- 
normalized airfoil shape to a full-sized airfoil as described 
above, then multiplied by a scaling factor equal to the 
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distance between the section locations. 



Section locations 



consisted of 24 wing-station positions along the leading-edge 
model, including the nine ribs and 15 intermediate points as 
determined by finite element model node locations. 

The first of the inputs to the program consisted of 
the airfoil geometry definition as described above. Next, a 
tabular file of wing station locations was read in together 
with. the section lift coefficient for that spanwise location 
as found by the span-load program. These section lift 
coefficients, initially found for the .35, .45 and .55t| 
locations (where t ]= 2y/b) were curve-fitted using the Cricket 

TM 

Graph plotting software in order to achieve a high degree of 
accuracy in determining the individual Cjg's at stations which 
were no more than three inches apart. Also included in this 
input file was a list of spanwise multiplication factors for 
scaling up the load as described above. Additional input 
files consisted of the finite element node numbers which were 
matched with their respective x and y direction loads in the 
program. Read in from the terminal were the airspeed under 
consideration and the twist angles of the wing box as 
determined from the span-load program. 

Output consisted primarily of the load file which included 
not only the loads at the finite element node points, but also 
the constraints and twist displacements for the upper flange 
and lower hinge node points of the leading-edge finite element 
model. The twist displacements were calculated within the 
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panel program using formulas based on geometric considerations 
of the height of the wing spar at the inboard and outboard 
ends of the leading-edge segment, and the net twist 
displacement of the front spar from the inboard to the 
outboard end of the model. That is, the finite element model 
was assumed to have undergone a rigid-body rotation to the 
degree of twist which was found to exist at the inboard end 
(WS 256). Therefore, twist displacements at the inboard end 
were set to zero, followed by application of the subsequent 
net twist distribution that occurred at the other wing 
stations while proceeding outboard to WS 320. This net twist 
distribution was equal to the difference between the outboard 
and inboard twist amounts, applied linearly over the 24 
stations. This twist amounted to approximately 0.3 degrees in 
the 275-knot, 3-G symmetric pullup. The displacements 
generated for application at the node points were in the 
longitudinal (x) and vertical (z) directions on the three- 
dimensional model. Twist rotation of the front spar was taken 
to be about its vertical mid point, and the structural wing 
box was assumed to have no chordwise distortion as it rotated 
about the .40c elastic axis location. Other output took the 
form of files to examine Cp distributions and to tabulate 
loads by wing station and two-dimensional node point for 
verification of the load file. In addition, output was 
generated which approximated the total normal and chordwise 
loads applied to the entire wing leading-edge center section 
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located between the engine nacelles. This estimation was 
accomplished by computing the load on the leading-edge portion 
of the airfoil at WS 274 and multiplying by the length of the 
leading-edge segment (92 inches). 

3 . Program Operation 

After reading in the coordinate information for the 
airfoil along with the wing station, Cjg and load 
multiplication factor, the FORTRAN program calculated the 
chord length at the particular station based upon the 0 . 40 
taper ratio. It also determined the thickness fraction of the 
airfoil section at that station by assuming a linear taper 
from 14% maximum thickness at the root to 12% at the wing tip. 
This thickness factor was then used to recalculate the y 
coordinate position of each node point to redefine the shape 
of the airfoil. An initial angle of attack of one degree was 
set and the process described in the theory section of this 
Chapter took place, wherein the Cy, and were calculated 
for that angle of attack. This Cjg value was then compared to 
that required (as input with the wing station), and an 
iterative cycle commenced in which the angle of attack was 
varied up or down by an amount based on the product of the 
deviation multiplied by a preset angular value. An accuracy 
test of .0001 was applied to reach an acceptable value for Cjg, 
at which time the process started over with the next wing 
station. Forces in the x and y direction were matched with 
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the appropriate finite element node points and written to the 
load file for each iterative cycle. After all loads were 
calculated and stored, the twist displacements and zero 
boundary constraints were calculated and appended to the load 
file. 

4. Verification of the Program 

The accuracy of the two-dimensional panel method was 
verified by comparison with published empirical data for 
tangential velocity and/or pressure coefficient distributions 
for the NACA 0012 and Eppler E64 airfoils before its use in 
this application. Results were nearly identical to the 
published data with only a small deviation seen near the 
trailing edge of the program tangential velocity distribution 
for the Eppler airfoil. No difference from the NACA 0012 Cp 
data could be identified. 

5. Flight Regime Selection 

It was found in the course of running the program at 
various airspeeds and angles of attack that the highest loads 
on the leading-edge segment were generated at slower airspeeds 
and higher angles of attack as a constant G load was 
maintained on the aircraft. This result seemed contrary to 
conventional opinion that the highest loads would most likely 
occur at or near the high speed end of the operating envelope. 
A brief study was undertaken to determine the cause of this 
phenomenon and a hypothesis is given here. 
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a. Method Employed 

In an effort to study the effect of dynamic 
pressure (airspeed) and angle of attack on the leading edge of 
an airfoil, the variation of net Cp distribution on the 
leading edge at various angles of attack was first examined. 
These Cp values were the numerical sum of the difference 
between the upper and lower Cp's, integrated over the 
chordwise distances occupied by their respective panels. 
These data obtained were then curve fitted with a third order 
polynomial and used in a spreadsheet to calculate the Cp 
distribution at varying angles of attack. These angles of 
attack were generated by varying the airspeed from 275 to 425 
knots, calculating dynamic pressure (q) for a 3-G wing loading 
from the Bernoulli equation, and then converting these q's to 
angles of attack required using the basic lift formula, 
altered by the equation. 



to give, 

a = SW/CL^qS and a = SW/CL^qS - C^o/C^a 
for the symmetrical NACA 0012 and cambered P-3 airfoil, 
respectively. These angles of attack were then used as inputs 
to the polynomial curve fits, from which a corresponding set 
of Cp values were calculated. Next the product of Cp and q 
were found, to give the pressure acting on the leading edge, 
corresponding to a matched set of q and angle of attack. This 
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information was then plotted as seen below in Figures 15 and 
16. 

b. Symmetric Airfoil 

Figure 15 shows the distribution of the various 
parameters on a relative scale as airspeed increases for the 
case of a symmetric airfoil . Note that the pressure acting on 
the leading edge is nearly constant, showing only a slight 
increase with increasing airspeed. 



Pressure Change of 0012 Airfoil LE with AOA and Q 




Figure 15. Pressure analysis of 0012 leading edge, 
c. Cambered P-3 Airfoil 

Next the same information is plotted for the P-3 
airfoil in Figure 16. Note that the pressure on the leading 
edge undergoes a marked decline as airspeed increases (AOA 
decreases) . 
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Pressure Change of P-3 Airfoil LE with AOA and Q 




Figure 16. Pressure analysis of P-3 leading edge, 
d. Analysis 

Because of the apparent drop in pressure on the 

cambered airfoil, it was concluded that the effect of camber 

was to shift the loading of the wing in a way that caused 

angle of attack to become the dominant influence rather than 

dynamic pressure. It was also observed that if the lift 

eguation applied in the program was that of a cambered airfoil 

(i.e. Cj 0 Q had some positive value as the lift curve was 

displaced upward), this loading phenomenon was present. If 

Cj 0 Q were zero, the load on the leading edge was approximately 

the same at various angles of attack and airspeeds. Because 

of this observation, it was decided to select the 275-knot, 3- 
* 

G flight position as the maximum loading position in the 
normal operating envelope . 
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V. 



FINITE ELEMENT ANALYSIS 



A. THEORY 

The method of finite element analysis is a means of 
simulating the structural behavior of a continuous physical 
system by a discretized representation of that system. 
Structures are represented by discrete node points which are 
connected by structural elements. The nodes form a grid which 
details the general shape of the structure while the elements, 
although they appear as only lines, are mathematically given 
the physical properties of the portion of the structure which 
they are there to represent. A physical structure is thus 
transformed into a mathematical representation for the purpose 
of analyzing some behavior of the structure. This analysis 
may be in the area of dynamic response, heat transfer, or, as 
is the case here, static loading response. This method of 
analysis is widely proven to be highly accurate and has been 
used in many engineering fields, including aerospace, 
automotive, civil and mechanical applications. With the 
increased capability, speed and data storage capacity of 
microcomputers, this analysis technique is no longer limited 
to mainframe applications, as was the case a few years ago. 
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1 . Accuracy 

The accepted rule of thumb in finite element modeling 
is that the use of more node points results in a more accurate 
solution. Convergence tables have been developed which show 
that the use of fewer nodes increases the stiffness of the 
model . The main drawback to using a large number of nodes is 
that it greatly increases computation time and requires larger 
amounts of storage space than a model of the same structure 
using fewer nodes [Ref. 14]. The finite element analysis 
software used for this thesis is called MSC/pal 2® and is a 
product of the MacNeal-Schwendler Corporation, the company 
which also creates the highly respected NASTRAN® finite 
element application for VAX/VMS work stations and mainframe 
operations. The accuracy of MSC/pal 2® has been tested and 
documented by the manufacturer [Ref. 15]. In addition, the 
manufacturer recommends simple hand calculations to verify 
that results obtained for a given model are reasonable (i.e. 
within the same order of magnitude). This was done by 
calculating simple beam bending stress with constant area 
cross sections approximating the rib legs of the model. 
Results established that the finite element model produced a 
solution which was well within expected norms. 
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2 . Equations 

The number of equations to be solved in the finite 
element analysis is equal to the number of degrees of freedom 
in the model. Each node point has six degrees of freedom: 
three translational (in each of the (x,y,z) coordinate 
directions, and three rotational (about each of the coordinate 
axes). Stiffness equations are generated for the stiffness of 
each connecting element, based on the specified material 
properties (Young's modulus, shear modulus, mass density, 
tensile yield stress are specified) and the geometric 
configuration of the element. Elements may take the form of 
beams, triangles, quadrilaterals and others. These nodal 
stiffnesses are combined to form a system stiffness matrix, 
[K], of size (NxN) where N is the number of equations. 
Degrees of freedom may be eliminated by setting them to zero 
in the model definition phase (a way of applying boundary 
constraints) or fully retained as was done in this 
application. (Here, boundary constraints were applied in the 
load as discussed in Chapter IV. This application allows 
variation of the boundary displacements from one load to the 
next and allows recovery of reaction forces at the constrained 
nodes.) Once the stiffness matrix has been formed, the static 
analysis may be performed according to the following equation: 

[K]{U} = {F} 
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where , 



• {F} = column vector of applied loads (Nxl) 

• {U} = resultant column vector of nodal displacements 

(Nxl) 

Gaussian elimination is employed to solve the matrix equation. 
In large models, as is the case here, matrix partitioning 
takes place prior to solution. 

Once the displacements are known, stress-strain 
relationships are employed to compute stress values throughout 
the structure. These stresses are available to the user in 
the form of major and minor principal axis stresses. Von Mises 
stress concentrations and maximum shear stresses. In 
addition, output of displacements and rotations are 
accessible. 

B. STRUCTURAL REPRESENTATION 

1. The Leading Edge Structure 

The leading edge segment considered for this analysis 
was the port wing, center section, located between the number 
three and four engine nacelles. Figure 17 shows a cutaway 
view of the structure. The segment is composed of 12 vertical 
ribs supporting a double (inner, outer) skin. The outer skin 
is .040 inches thick and the inner is a stamped corrugation of 
.016 inches. Assembly of the outer and inner skins provides 
a series of ducts approximately .25 inches in height and two 



49 




Figure 17. P-3 leading edge structure [Ref 16]. 

inches wide each. These ducts run longitudinally along the 
inner surface of the leading edge and are continuous over the 
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spanwise length of the leading edge. The purpose of this 
series of ducts is to provide a channel through which bleed 
air may travel to heat the leading edge. The tube seen 
extending through the structure delivers the bleed air to the 
ducts. The leading edge segment is secured to the front spar 
of the wing by a full-length piano hinge at the bottom edge 
and a screwed-down spar cap flange at the top. This 
arrangement provides access to the area for maintenance 
functions . 

2 . The Model 

The basic finite element model used in this analysis 
was obtained through translation of a NASTRAN® finite element 
model using a function in the MSC/pal 2® application called 
NASPAL®. NASPAL® reads the NASTRAN® text file for the model 
and rewrites it in the format used by MSC/pal 2®. The 
NASTRAN® model file was obtained from Aerostructures, Inc. 
through contact with NAVAIR's AIR-530 office. Because the 
NASTRAN® model consisted of 2749 nodes, it was necessary to 
reduce the model size to meet the MSC/pal 2® limitation of 
2000 nodes. This reduction was done by entering a set of 
geometric coordinates during the NASPAL® translation and 
instructing the translator to consider only the outboard 64 
inches of the model. In effect, the leading edge section was 
severed between WS 247 and WS 256, or between the third and 
fourth ribs from the left end as shown in Figure 17. The 
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remaining portion of the model was translated from the fully 
defined model. This approach was considered to provide a more 
accurate solution than increasing the spacing between nodes, 
keeping the highest available 
level of detail in the model 
(again, more nodes mean 
better accuracy for the same 
structure). In doing so, it 
was necessary to modify the 
inboard end rib of the 
structure (WS 256) since the 
end ribs were constructed 
differently than the 
intermediate ribs. The end 
ribs are closed in the front 
and have single, riveted 
flanges on their interior 
(Figure 18, bottom), while 
the intermediate ribs (Figure 
18, top) have an open front 
to allow access for the bleed 
air "pump cap" assembly which 
leads to the double skin 
described earlier. The 

intermediate ribs also have 
double flanges as shown. 




Figure 18. Intermediate and 
end rib detail. 
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These changes to the model were made in the rib at WS 256 (now 
the inboard end rib of the model) by adding node points and 
quadrilaterals to close the front and by rearranging the beam 
element properties within the MSC/pal® 2 model definition text 
file. In this way, the model was altered to represent a 
shortened leading edge segment with properly defined ribs. 
The model was constructed as an all-aluminum structure, and as 
stated in Chapter IV, the loads applied were generated for 
each specific wing station and scaled according to the 
spanwise distance between the node points. In this way, the 
model detailed here received a scaled-down load for its 
scaled-down size. The upper and lower surface views of the 
model are presented in Figures 19 and 20. The end views of 
the model are shown in Figure 21. 

The lower hinge is replicated in the model by leaving 
the Y-axis rotation unrestrained. The upper flange of the 
model is secured in all six degrees of freedom. Displacements 
for front spar twist are incorporated in the upper and lower 
constraints as detailed in the previous Chapter. All of these 
boundary conditions are input through the "Displacements 
Applied" command section of the load file. 

Another difference between this model and the original 
NASTRAN® code developed for NAVAIR is that the original did 
not incorporate stiffness generation in the skin of the model. 
The skin thickness of the NASTRAN® model was .056 inches, 
which is the sum of the outer and inner skin thicknesses but 
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Figure 19. Upper surface of finite element model. 
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without the 
stiffness 
generation enabled 
in the model 
definition file, 
the skin would have 
no stiffness. In 
effect, it would 
act as a non- 
loadbearing 
membrane. In the 
MSC/pal 2® model, 
stiffness 
generation was 
enabled, but when 
loaded, this 
resulted in 
excessive 
deformation of the 
skin. It was 

decided that the 
.056 inch thick 
skin did not accurately model the combined effect of the inner 
and outer skin combination since the corrugated inner skin 
would be much stiff er than its mere thickness (.016 inches) 
would represent. Although there is a variable stiffness 




Figure 21. End views of finite element 
model (outboard, inboard). 
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factor in the MSC/pal 2® application, available references on 
finite element analysis did not detail the calculation of this 
factor. Therefore, it was decided to increase the skin 
thickness to an equivalent thickness which would accurately 
replicate the behavior of the corrugated skin combination. 
This thickness was calculated by equating the moment of 
inertia of a box beam formed by a single corrugation and outer 
skin combination with that of a solid cross section about the 
same reference axis (at the upper surface). After allowing 
for spaces between the corrugations where the two skins are 
riveted together, an equivalent skin thickness of .1838 inches 
was determined and employed. This action reduced the 
deformation of the skin under load to reasonable norms. 

Loads were applied to the model at spanwise rows of 
selected node points which most nearly corresponded to the mid 
points of the panels in the two-dimensional panel method. 
These point loads were applied in the vertical (Z) and 
horizontal (X) directions in units of pounds force. 

Material properties used to represent the 2024-T6 
aluminum structure in the construction of the model were as 
follows: 

• Young's modulus (E) = 1.06E+07 psi [Ref. 17] 

• Shear modulus (G) = 4.0E+06 psi [Ref. 17] 

• Poisson's ratio (u)= 3.25E-01 

• Tensile yield strength (c7yg) = 4.7E+04 psi [Ref. 18] 
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where Poisson's ratio (u) was calculated according to the 
standard formula: 



G= 



E 

2(l+v) 



C. APPLICATION 

Many load conditions were examined in the course of this 
thesis. Presented here are the six load cases which give the 
best overall illustration of the observed effects of static 
aerodynamic loading, both with and without static aeroelastic 
effects included. In this way, the effect of wing box twist 
may be seen, along with the combined effect of angle of attack 
and dynamic pressure. In Table 3, the features of the six 
load cases are given. The L/R column provides a distinction 
between a wings level (L) pullup or a rolling (R) pullup. The 
effect of rolling into a turn during the application of G 
loading (as in a climbing breakaway maneuver) is not the same 
as that of rolling out of a turn during G application (as in 
rolling to wings level while pulling out of a dive) . Since it 
was found that the loading effect in terms of both twist and 
air loading was greater during the former (due to the combined 
effect of roll helix angle, aileron deflection and air load), 
the former was chosen for presentation here in the rolling 
load cases. All load cases were generated at 135,000 pounds 
gross weight except for number 4 which was done at 110,000 
pounds. In load case number 1, the effect of twist was 
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eliminated from the load solution by setting the MSC/pal 2® 
displacements to zero so that the effect of wing torsional 
twist may be seen by comparison with load case 2. 



TABLE 3. LOAD CASES EMPLOYED 



Load Case 


Airsoeed^ kts^ 


Load Factor 


L or R 


Comment 


1 


275 


3.0 


Level 


No twist. 


2 


275 


3.0 


Level 




3 


275 


2.4 


Roll 




4 


240 


3.0 


Level 


110,000 # 


5 


350 


2.4 


Roll 




6 


325 


4.5 


Level 





1 . Input Data 

Presented in Table 4 are the input data sets for each 
of the load cases. These are given for the three pertinent 
non-dimensional spanwise wing stations (2y/b) as generated by 
the span load program. The 24 actual wing stations (2y/b = 
0.4322 through 0.5397) employed for the two-dimensional panel 
program were solved by curve fitting the lift coefficients 
bounded by these extremes, as described in Chapter IV. The 
92-inch loads are the approximate longitudinal (positive aft) 
and vertical (positive up) total loads which would be seen by 
a complete center section leading edge segment on the 
aircraft. The total load applied to the shortened finite 
element model would be some 69.6 percent of this stated load. 
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TABLE 4. FINITE ELEMENT ANALYSIS INPUT DATA 



Load Case 


2y/b 


Ci 


Twist ^dea) 


92" Load (lbs 


1 


.35 


1.2726 


0.0000 


FX = -5157 




.45 


1.2379 


0.0000 


FZ = 13155 




.55 


1.1834 


0.0000 




2 


.35 


1.2726 


0.8736 


FX = -5157 




.45 


1.2379 


1.1221 


FZ = 13155 




.55 


1.1834 


1.3764 




3 


.35 


1.0778 


0.7605 


FX = -3657 




.45 


1.0569 


0.9966 


FZ = 11069 




.55 


1.0090 


1.2572 




4 


.35 


1.3465 


0.7216 


FX = -4454 




.45 


1.3118 


0.9272 


FZ = 10666 




.55 


1.2575 


1.1375 




5 


.35 


0.7109 


0.5673 


FX = -2103 




.45 


0.6723 


0.7335 


FZ = 10540 




.55 


0.6100 


0.9169 




6 


.35 


1.3789 


1.2997 


FX = -5982 




.45 


1.3247 


1.6613 


FZ = 16742 




.55 


1.2453 


2.0279 





2. Finite Element Analysis Results 
a. Load Case 1 

In looking at the 275 knot load without static 
aeroelastic twist, it was noted that the largest stress 
concentrations in the structure were located in the rib legs, 
with the lower legs experiencing approximately 15 ksi in 
tension along their upper flanges and the upper legs seeing 
cibout -14 ksi (compression) along the lower flanges. These 
stresses were evenly distributed, as may be seen by the values 
in Figures 22, 23 and 24 where the major principal axis stress 
(Oj) contours are shown. The minor principal (ajj) stresses. 
Von Mises and shear stresses show a similarly even 
distribution, with all stress levels well below the yield 
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stress value (o^g = 47 ksi) for the material. It should be 
noted that the apparent deformations in the plot are 
exaggerated and not to scale. This scaling provides the 
viewer with a better perception of the direction of 
displacement occurring, although in reality, the displacements 
are only on the order of 0.06 inches for this load case as 
determined by MSC/pal 2®. 










MAJOR STRESS 

A-2.2483E*03 
B 0.0000E^00 
C 2.2483E+03 
D 4.4966E+03 
E 6.7449E+03 
r 8.9932E+03 
G 1.1241E+04 
H 1.3490E+04 
! 1.5738E+04 
J 1.7986E+04 



Figure 22. Major principal axis stress contours on 
inboard (WS256) end rib in untwisted 
condition (Case 1). 
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B 0.0000E400 
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E 6.7449E+03 
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I 1.5738E+04 
J 1.7986E+04 



Figure 23. Major principal axis stress contours on 

middle (WS282) rib in untwisted condition 
(Case 1). 
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B 0.0000E«00 
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H 1.3490E+01 
I 1.5738E+01 
J 1.7986E+04 



Figure 24. Major principal axis stress contours on 



outboard (WS320) rib in untwisted condition 
(Case 1 ) . 
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The untwisted load case (case 1) does not 
accurately reflect the behavior of the leading edge since it 
does not include the effect of aeroelastic-induced spanwise 
twisting (torsion) of the wing. The stress contours in 
Figures 22, 23 and 24 were presented for comparison to load 
case 2. 

i>. Load Case 2 

A dramatic difference in the observed stress 
contours occurred when the finite element model was subjected 
to the net twist occurring in the wing box as determined by 
the static aeroelastic span load analysis. It should be noted 
that the displacement input to the model equated to only about 
0.3 degrees of front spar rotation from WS 256 to WS 320. The 
same set of stress contour plots as above are given in Figures 
25, 26 and 27. A striking contrast existed between the first 
and second load cases, with the case 2 major principal axis 
stress distribution very unevenly spread between the inboard 
and outboard ends. As may be seen in Figure 25, the inboard 
end rib (WS 256) experienced more than double the stress 
concentration in its lower leg with 34.8 ksi being the highest 
level contour shown. Moving outboard, WS 282 (Figure 26) 
showed a maximum of about 20 ksi in its lower leg while the 
stress in the same leg on the outboard (WS 320) rib went to 
zero. (A stress level of 2.9 ksi in the saddle and upper leg 
of the rib may still be seen. ) 
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H 3.4842E+04 
0 3.7745E+04 



Figure 25. Major principal axis stress contours on 

inboard (WS256) end rib with twist applied 
(Case 2). 
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Figure 26. Major principal *axis stress contours on 
middle (WS282) rib with twist applied 



(Case 2 ) . 



63 



MAJOR STRESS 



Z 




A-2.9035E+03 
B 0.0000E«00 
C 2.9035E+03 
D 5.8069E«03 
E 8.7104E+03 
r 1.1614E+04 
G 1.4517E+04 
H 1.7421E+04 
I 2.0324E+04 
J 2.3228E+04 
K 2.6131E+04 
L 2.9035E+04 
n 3.1938E+04 
H 3.4842E+04 
0 3.7745E+04 



Figure 27. Major principal axis stress contours on 
outboard (WS320) rib with twist applied 
(Case 2). 



Since it was observed that this same pattern of 
stress distribution occurred for the minor principal axis, Von 
Mises and maximum shear stresses, only the inboard end (WS 
256) contours are shown in Figures 28, 29 and 30. This 
pattern is to be expected since they are all geometrically 
related. The major and minor principal stresses occur on 
planes on which there is no shear stress and are oriented 
perpendicularly to each other, while the maximum shear stress 
occurs on planes which are at angles of 45 degrees to the 
principal planes. Von Mises stresses (a^) are derived from a 
criterion known as the Maximum Distortion Energy Criterion and 
may be found from the equation, 

V = 
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which is based on the determination of the energy associated 
with changes in shape of a given material. These relation- 
ships are valid under the assumption of a plane stress 
condition in the material. This means that the metal is thin 
in comparison to its other dimensions so that the stresses 
across the thickness of the metal may be considered as 
negligible. This plane stress condition is the case for most 
aircraft structural components since they are made of sheet 
material, and is valid here. [Ref. 19 , 20 ] 
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Figure 28. Minor principal axis stress contours on 

inboard (WS256) end rib with twist applied 
(Case 2 ) . 
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Figure 29. Von Mises Criterion stress contours on 

inboard (WS256) end rib with twist applied 



(Case 2 ) . 




M8X SHR/ZND NJ 

ft 0.0000E«00 
B 1.8091E«03 
C 3.6183E*03 
D 5.4274E*03 
E 7.2366E*03 
r 9.0457E*03 
G 1.0855E«04 
H 1.Z66^E404 
I 1.4473E*01 
J 1.6Z8ZE«04 



Figure 30. Maximum shear stress contours on inboard 

(WS256) end rib with twist applied (Case 2). 



In Figures 31 through 36, the x and z displacements of 
the ribs are shown in pairs from inboard to outboard. Note 
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the small displacements occurring in the upper flange and 
lower hinge due to twist of the front spar. 
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Figure 31. Displacement (x) at inboard rib (WS256). 
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Figure 32. Displacement (z) at inboard rib (WS256). 



67 



ccc 




TRANS. DEFL. X 

fl-3.6889E-02 
B-2.7667E-02 
C-1.84i4E-02 
D-9.2222E-03 
E 0.0000E400 
r 9.2222E-83 
G 1.844iE-82 
H 2.V667E-82 
1 3.6889E-82 
J 4.6111E-82 



Figure 33. Displacement (x) at middle rib (WS282). 
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Figure 34. Displacement (z) at middle rib (WS282). 
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Figure 35. Displacement (x) at outboard rib (WS320). 
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Figure 36. Displacement (z) at outboard rib (WS320). 



Another available representation of the stress 
concentrations present in the structure is an X-Y plot of all 
four stresses together in a form resembling a frequency 
scatter on an oscilloscope. This plot for the case 2 load 
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condition is presented in Figure 37. While this is a rather 
cluttered plot it does serve to provide, at a glance, the 
maximum and minimum stress concentrations in the structure. 
Since the stress contours of individual members in the 
structure have been examined and it has been found that the 
highest stresses are in the rib legs, the peak stresses 
depicted may be attributed to these locations. The positive 
(upper) portion of the plot is a combination of the major 
principal. Von Mises and shear stresses while the lower half 
depicts the minor principal axis stresses. For the sake of 
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Figure 37. X-Y plot of combined stresses (Case 2). 



brevity, and since the contour plots of the individual ribs 
retain the same relative form as those presented for cases 1 
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and 2, the stress levels for the remaining load cases will be 
presented in this form. 

c. Load Case 3 

For the 275-knot, 2.4-G rolling pullup, the 

maximum stress levels observed are approximately 33 ksi on the 
descending wing as shown in Figure 38. The loads on the 
leading edge of the ascending wing are lower than those for 
the descending wing. This speed was chosen because of the 
relationship which was shown to exist between leading edge 
loading, angle of attack and dynamic pressure. The result of 
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Figure 38. X-Y plot of combined stresses (Case 3). 
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higher airspeed at the same aircraft load factor may be see in 
part e of this Chapter. 

d. Load Case 4 

The 110,000 pound weight condition was chosen for 
this load case in order to examine the reduction in leading 
edge stress when operating within the normal flight envelope 
at less than maximum gross weight. The 240-knot airspeed is 
the approximate "corner" speed at a load factor of 3G for this 
weight. The maximum stress values seen under this condition 
are 28.8 ksi in tension and 28.3 ksi in compression as seen in 
Figure 39. 
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Figure 39. X-Y plot of combined stresses (Case 4). 
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e. Load Case 5 



In load case 5, the effect of a 2.4-G rolling 
pullup at 135/000 pounds gross weight at a true airspeed of 
350 knots was examined. In Figure 40, it may be seen that the 
maximum stress values are approximately 26 ksi, which is 
approximately 7 ksi less than the same maneuver at 275 knots 
( load case 3 ) . 



350 KT. 2HG. ROLLING PULLUP 
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Figure 40. X-Y plot of combined stresses (Case 5). 
f. Load Case 6 

In order to examine the design ultimate flight 
load condition, load case 6 was run at 325 knots and 4.5G. 
The airspeed corresponds to the approximate stall speed 
(corner speed) for the 135,000 pound airplane gross weight. 
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Stress values of nearly 52 ksi may be seen in Figure 41. Note 
that this result exceeds the yield stress for the material (47 
ksi) . 
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Figure 41. X-Y plot of combined stresses (Case 6). 



3. Discussion 

The primary factor of concern in the observed loads . 
presented above is the apparently high level of stress in the 
rib legs as determined by the methods stated in the course of 
this writing. According to the Engineering Investigation (El) 
report for the USN P-3 mishap [Ref. 21] and a preliminary 
report issued by Aircraft Research Laboratories of Melbourne, 
Australia [Ref. 22], the initial structural failure in both 
mishaps was believed to have occurred in the lower rib legs. 
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with initial failure occurring in the outboard leg and 
proceeding inboard in a series of sequential failures. It 
would appear from the results of this study that this failure 
sequence may have been reversed, since the highest observed 
stresses were in the inboard rib. Additionally, the fact that 
there is an intermediate rib at WS 317, only three inches from 
the outboard end rib (rather than the eight- to nine-inch 
spacing for all other ribs in the leading edge) supports the 
idea that the failure would likely have initiated at some 
other location. This close proximity of two ribs means that 
the load in that region would be shared, thereby reducing the 
stress in the outboard rib. In any case, the primary question 
remains as to whether there were sufficiently high stress 
levels in the rib legs to cause structural failure. 

While the maximum observed stress levels (found here 
for operation within the flight envelope, occurring under load 
case 2 in a wings level, 3-G pullup at 135,000 pounds gross 
weight) were some 12 ksi below the yield stress for the 
material, two areas of concern must be addressed: 1) the 
additional effect of extending this data from the assumed 64- 
inch model size to the true 92-inch leading edge segment which 
is installed on the aircraft; and 2) the effect of stress 
concentrations around holes in the rib legs. 
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a. Extending the Structure 

Proceeding under the assumption that the stress 
concentrations vary linearly with the length of the structure, 
extending the 64-inch model to the full (92-inch) length of 
the true structure would result in an increase in the observed 
stress in the lower rib leg to 51.9 ksi, which is in excess of 
the 47 ksi yield stress. Examining the major principal axis 
stresses in the lower legs of WS 256 and WS 282, this 
linearity would appear to hold true. Since aluminum is a 
ductile material, the stresses may subside as the material 
begins to yield, therefore, the structure may still carry the 
load without fracture. However, changes in the pressure 
distribution around the leading edge must also be considered 
as the structure begins to deform. Given the observed 
directions of the deformations as seen in the contour plots, 
it seems likely that the aerodynamic load would continue to 
increase due to the increased camber in the leading edge as it 
lifts perpendicular to the chord and develops a more circular 
form. This deformation, although initially small could add to 
the load, initiating a divergent stress scenario leading to 
failure. 

Jb. Stress Concentrations 

It is an established fact that stress 
concentrations around circular holes can multiply the 
localized stress by factors of between two and four, depending 
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on the directional orientation of the loading seen by the 
specimen. Along the inner surface of the end ribs is a 
continuous flange which is attached by a rivet spacing of 
approximately one-half inch. Additionally, there are four 
tapered angle stiffeners riveted to the web section of the 
rib. Two of these are symmetrically located on the upper and 
lower legs at a position approximately 52 percent of the 
leading edge chordwise dimension forward of the hinge. (On 
the outboard end rib (WS 320) this is approximately 12 inches 
forward of the hinge.) Around these rivet holes, stress 
concentrations may be assumed to occur. The occurrence of 
these concentrations at the onset of stress loading depends on 
the rivet installation procedure; i.e., the stress 
concentration effect is delayed if the rivet hole is placed in 
compression (as in a dimpled or double dimpled installation 
[Ref. 18]). If the hole were initially in compression, a 
margin or stress buffer would be provided, as the material is 
subjected to tension stress; i.e., the tensile stress applied 
must first exceed the compressive pre-stress of the rivet 
installation before the hole will begin to experience a stress 
concentration. The type of installation of the rivets under 
consideration is unknown by this author. Since the inside 
flange of the upper legs is placed under compressive loading 
nearly equivalent to that experienced by the lower leg in 
tension, the effect of compressive concentrations must also be 
considered. If the rivet holes were pre-stressed in 



77 



compression, this would add to the compressive stress level at 
this location, potentially causing the upper leg to fail first 
in compression. As before, this all depends on the stress 
reduction due to yielding which takes place in the ductile 
material . 

In the Australian incident report, the failure 
sequence was described as having initiated along the line of 
the riveted stiffener on the lower leg of the outboard end rib 
as mentioned above. In addition, it goes on to say that the 
failure then proceeded along a line of dimples in the lower 
legs of the intermediate ribs, from outboard to inboard. 
These dimples act as stiffeners for the intermediate ribs and 
would also act as stress concentrations under load. None of 
these stiffeners, on either outboard or intermediate ribs, nor 
their rivet holes, are included in the finite element model 
employed in this study. Again, the failure sequence described 
is opposite of that expected from the results of this study. 
The RAAF determination of the failure sequence was based on 
the observed outboard-to-inboard bending of the rib fragments 
which remained attached to the structure. [Ref. 22] 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

The results of this study indicate that failure of a P-3 
wing leading edge segment could be predicted to occur within 
the normal aircraft operating envelope providing proper 
account is given to the influences of wing static aeroelastic 
effects both upon wing span loads and torsional twists induced 
in the wing spar box. The location of the highest stresses 
are in the area of the observed failures as reported in the El 
and Preliminary reports [Ref. 21 and 22]. It is reasonable to 
assume that stress concentrations around the rivet holes could 
be on the order of 1.5 or 2.0 times the observed stresses. 
This, combined with the effect of extending the model to the 
full 92-inch length, could result in stress levels in excess 
of the ultimate strength of the material, 60 ksi [Ref. 18]. 

The primary evidentiary conflict with this study is the 
reported location of apparent failure initiation at the 
outboard vice inboard end. The observed inboard bending of 
the rib fragments may have been induced by the upper portion 
departing the wing in some fashion other than that assumed. 
It is difficult, at best, to assess the direction in which the 
ribs would bend as the failure progressed. 
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While fatigue is potentially a contributor to this mode of 
failure, given the possibility of plastic strain having 
occurred earlier under a reduced load condition, the El [Ref. 
21] stated that no undue hardness was detected in the failed 
USN structure. This would indicate that no material strain 
hardening had taken place up to the time of failure. 

Corrosion is another potential problem, particularly when 
dealing with aircraft operating routinely in low level 
overwater environments; however, no levels of corrosion which 
would effect the structural integrity of the leading edge were 
cited in the El. 

Whether or not the aircrews overstressed the aircraft by 
exceeding the limitations of the flight envelope cannot be 
concluded here. Even they may not know for sure whether this 
was the case. Flight experiments conducted in the 2F87F P-3 
simulator at NAS Moffett Field, CA, during the course of this 
analysis indicate that at high gross weights and aft center of 
gravity conditions, it is easy to exceed the 3G limit with 
just a firm pull on the control yoke. The rapid onset of G 
overload may be imperceptible due to the rate at which it was 
observed to occur. The cockpit G meter is not an instrument 
which is normally kept in the pilot's scan. 



80 



B. RECOMMENDATIONS 



This analysis should be confirmed by independent 
validation. If these results are verified, re-enforcement of 
the end rib may be feasible; alternately, some limitations 
should be placed on the operational flight envelope of the 
P-3 . Covered here are some recommended procedures for 

validation along with approximations of interim flight 
envelope limits, should validation and/or repair not be 
possible in the near-term. 

1. Validation 

Since the Naval Air Systems Command, through 
Aerostructures, Inc., has the base model from which this 
finite element model was derived, it is recommended that 
validation be conducted with the following modifications on 
the existing model: 



• Skin stiffness should be included in the model. The 
contribution of the skin to the load bearing properties of 
the structure is essential. Skin stiffness could be 
simulated by artificially increasing the thickness to 
simulate the double skin as was done here, or preferably, 
through addition of a w axis (quadrilateral element local 
coordinate system) stiffness factor which would emulate 
the combined stiffness of the inner and outer skins. 

• Static aeroelastic twist must be included in the analysis 
since it was the factor contributing to the dramatic 
increase in inboard end rib stress. Since the NASTRAN® 
model is already of the appropriate length, this model 
would serve to verify or discount the theory that the 
stress concentrations would increase linearly with 
extended length. 

• A small finite element model of the area of immediate 
concern (the lower leg of an end rib) should be 
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constructed and analyzed. This model should be of 
sufficient detail to include the rivet holes and 
stiffeners in order to examine the amount of stress 
concentration occurring. Any pre-stress of the rivet 
holes should be taken into account. 



Laboratory tests of actual end ribs would prove 
helpful in determining the effect of the rivet holes on stress 
concentrations and provide data on the load levels required to 
cause fracture. 

2. Flight Envelope Modifications 

If these results are verified, the P-3 flight envelope 
should be restricted to prevent another mishap. If, in the 
estimation of NAVAIR, verification cannot be accomplished 
within a reasonable amount of time (six months), then it is 
recommended that an interim measure be taken to restrict the 
flight envelope until such verification or disproval can be 
accomplished. While the precise envelope restrictions would 
be developed by NAVAIR, it is recommended on the basis of the 
results observed here that the following limitations be 
applied if interim limitations are deemed appropriate: 

• Reduce the maximum sustained load factor for the aircraft 
from 3G to 2G at all operating weights above 100,000 
pounds . 

• Reduce the maximum sustained load factor for the aircraft 
from 3G to 2.5G at all operating weights up to and 
including 100,000 pounds. 
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While these are rough estimates, they are believed to provide 
a sufficient margin of safety to allow safe flight without 
severely impacting daily P-3 operations. 
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APPENDIX A 



10 ' Program »‘P3SPNLD3 . BAS" 

' Solve Static Aeroelastic Spanload Problem 
' for a wing with straight taper 

' Vortex Lattice capability uses swept bound vortex 
' elements 
' Inputs consist of: 

' AR = Wing Aspect Ratio = B^2/S 

' TR = Wing Taper Ratio = Ct/Cr 

' SWP25 = 0.25 ch'd. sweep angle, +'ve is sweepback 
' CEA = Elastic axis location on const, fraction 
' chord line 

' MACH = Subsonic Mach no. for aerodynamic compress. 

' correctn. 

' M = No. of equal length spanwise stations, RH wing 

' N = No. of equal chordwise stations 

' ** Comments ** 

' Allows Vortex Lattice solution with MxN boxes on 
' RH wing 

' Max. of MxN =50:' Consistent with dimension 
' statements 

' N = 1 for elementary (Modified Weissinger) lifting 
' line theory 

' Wingspan, B, set to 1188 inches (for P3 applicatn. ) 
100 ' $DYNAMIC 

DIM Xl(50), Yl(50), X2(50), Y2(50), X3(50), Y3(50), 
SWP(50),DIM Al(50, 50), A2(50, 50), ASYM(50, 50), 
AANT(50, 50), DIM S(50, 50), F(20, 50), G(20, 50), 
XEA(20), YEA(20), El (20), GJ(20),DIM SWPM(50), 

DELA(50, 20), MX(20, 50), MY(20, 50 ) , DIM ALPHA( 50 ) , 

U(20, 20), UGJ(20, 20), UEI(20, 20),DIMA(50, 50), 
XIN(50): ' for Spanload solutns. using ELU/SLVB S/R's . 
DIM ALPHA1(50), ALPHA2(50), ALPHAS (50), ALPHA4(50), 
Q(50), DELCAM(5) 

OPEN "C: QBFILES\EXTRA.DAT" FOR OUTPUT AS #1 
150 ' Input wing geometric information . . P3 usage 
AR = 7.539: TR = .40088: M = 10 
'PRINT : PRINT "Aspect Ratio, AR ="; AR 
'PRINT "Taper Ratio, TR ="; TR 
'PRINT "No. R.H. Wing Spanwise Stas. ="; M 
'MACH = .6: PRINT TAB (5); "Mach No. ="; MACH 
INPUT "Airspeed (knots) ="; VEL 
MACH = (VEL * 1.68781) / 1116.3 
PRINT "Mach No. =", MACH 
IF MACH > .95 THEN GOTO 151 
FC = 1! / SQR(1! - MACH " 2): GOTO 152 
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151 FC = 11 

152 ’ Continue Dummy statement 
N = 5 

'PRINT "No. Chordwise Stas. ="/ N 
'INPUT "No. Chordwise Stas. ="; N 

SWP25 = -1.312: 'PRINT "C/4 Sweep (deg.) ="; SWP25 

'INPUT "0.25 chd. sweep (deg.) ="; SWP25 

CEA = .4: 'PRINT "Elast. Axis Locatn. , X/C ="; CEA 

'INPUT "Elastic Axis Locatn. , X/C ="; CEA 

'INPUT "Dynamic Pressure, Q (psi) ="; Q 

INPUT "Is this a RIGID run? l=Yes, 2=No, ANS=", ANS 

IF (ANS = 1) THEN 

Q = 0 

ELSE 

Q = (.7 * 2116.2 * MACH " 2) / 144 
END IF 

PRINT "Q =", Q 

KMAX = M * N 
/ 

' Calculate propeller thrust per engine (using curve fit 
' from 2F87F P-3 simulator data) 

THRUST = (5849 - 5.069 * VEL) 

' Calculate velocity boost due to prop slipstream 

V = VEL * 1.68781 

VI = (-V + SQR((v '' 2 + (4 * .05 * THRUST / (143.14 * 
.0023769))))) * 10 

' Use velocity boost to find boosted Q (Ql) 

Q1 = (.5 * .0023769 * (v + VI) " 2) / 144 
IF (ANS = 1) THEN 
Ql = 0 
END IF 

' Set up Q as a vector with boosted pressure (Ql) in 

' region behind propellers 

FOR I = 1 TO KMAX 

Q(I) = Q: NEXT I 

FOR I = 11 TO 35 

Q(I) = Ql: NEXT I 
/ 

PI = 41 * ATN(1!): ' Establish constant 

SWP25 = SWP25 * PI / 1801: ' Convert sweep to radians 

180 ' Print header 
'GOSUB 1100 
'FOR MM = 1 TO 4 



200 ' Determine Initial Wing Geometry 

B = 1188!: ' DEFAULT Value of P3 Wing Span, inches 
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DELB = .5 * B / M; 'Spanwise spacing increments of 
' vortex elements 

CR = 2 ! * B / (AR * (11 + TR) ) : ' Root Chord, inches 
S = B ^ 2 / AR: ' Wing Area, sq. in. 

TANLE = TAN(SWP25) + (.5 * CR * (11 - TR) ) / B: 

' Tangent L.E. sweep 

' Develop Mean Aero. Chord information 

MAC =2! * CR * ((11 + TR) ~ (TR / (11 + TR) ) ) / 31: 

' Cmac , inches 

IF TR = 11 THEN GOTO 210 

YMAC = .5 * B * (11 - (MAC / CR) ) / (11 - TR) : GOTO 211 

210 YMAC = .5 * B 

211 XMAC25 = (YMAC * TANLE) + .25 * MAC 

250 ' Determine coords, for wing vortex lattice corners & 

' control pts. 

CONSTl = 11 - TR 
FOR I = 1 TO M 

Cl = CR * (11 - (CONSTl * (I - 11) / M)): 'Total wing 

' chord, inbd. vortex 

C2 = CR * (1! - (CONSTl * I / M)): ' Wing chord, outbd. 

' vortex 

C3 = CR * (11 - (CONSTl * (I-.5)/M)): ' Chord at 

' control point sta. 

DELCl = Cl / N: DELC2 = C2 / N: DELC3 = C3 / N 
YEA(I) = (I - .5) * DELB: ' Create "M” Elastic axis 

' coordinates 

XEA(I) = YEA(I) * TANLE + CEA * C3 
FOR J = 1 TO N 

K=(I-1)*N+J: 'Create Vortex Lattice numbering 
' scheme 

Y1(K) = (I - 1) * DELB: Y2(K) = Y1(K) + DELB 
Y3(K) = Y1(K) + (.5 * DELB) 

X1(K) = DELCl * (J - .75) + Y1(K) * TANLE 

X2(K) = DELC2 * (J - .75) + Y2(K) * TANLE 

X3(K) = DELC3 (J - .25) + Y3(K) * TANLE 

TANSWP = (X2(K) - X1(K)) / (Y2(K) - Y1(K)) 

SWP(K) = ATN( TANSWP) 

TANSWPM = FC * TANSWP 
SWPM(K) = ATN( TANSWPM) 

NEXT J 
NEXT I 

' Tangent of Elastic Axis Sweep, etc. 

DELEA = (XEA(2) - XEA(l)): REA = SQR(DELEA " 2 + DELB"2) 
TANEA = DELEA / DELB: SINEA = DELEA / REA 
COSEA=DELB/REA 
WRITE #1, Q 

INPUT "ENTER G LOAD, n =", NZ 
'Calculate required CL 

CLREQD = (NZ * 135000) / (.7 * 2116.2 * MACH " 2 * 1300) 
PRINT "CLREQD =", CLREQD 
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270 



WRITE #1, NZ 

' Bring in El (I) and G J ( I ) values for Structural model 
'1=1 (Root Sta.) & I = M (Tip Sta.)f M = 10 by 
' default 



' El Data, Est 
DATA 8.30E+10, 
DATA 2.00E+10, 
FOR I = 1 TO M 
' GJ Data, Est 
DATA 7.50E+10, 
DATA 1.30E+10, 
FOR I = 1 TO M 



for P3 wing 
6.80E+10, 5.27E+10, 
1.33E+10, 0.85E+10, 
READ El (I): NEXT I 
for P3 wing 
5.25E+10, 3.90E+10, 
0.85E+10, 0.65E+10, 
READ G J ( I ) : NEXT I 



4. lOE+10, 
0.55E+10, 



2.90E+10, 

0.35E+10, 



2.80E+10 

0.35E+10 



2. lOE+10 
0.25E+10 



300 'Develop Aerodynamic Influence Coefficients 
' A1(K1,K2) = RH Wing 
' A2(K1,K2) = LH Wing 

FOR K1 = 1 TO KMAX: ' K1 Sta. is at the Control Point 

FOR K2 = 1 TO KMAX: ' K2 Sta. is at the Vortex Station 

NUMl = X3(K1) - XI (K2): NUM2 = X3(K1) - X2(K2) 

NUMl = FC * NUMl: NUM2 = FC * NUM2 

' Bring in Prandtl-Glauert factor 

DENI = Y3(K1) - Y1(K2): DEN2 = Y3(K1) - Y2(K2) 

R1 = SQR(NUM1 " 2 + DENI " 2 ) : R2 = SQR(NUM2"2+DEN2"2 ) 
' Find trig, functns for orthogonal transformation on 
' swept bound vortex 

SINSWP = SIN(SWPM(K2) ) : COSSWP = COS ( SWPM(K2 ) ) 

H = NUMl * COSSWP - DENI * SINSWP 

YIROT = NUMl * SINSWP + DENI * COSSWP 

Y2ROT = NUM2 * SINSWP + DEN2 * COSSWP 

COSTHETl = YIROT / Rl: COSTHET2 = -Y2ROT / R2 

' Logic Check to avoid division by zero 

IF (ABS(H)) <= .001 THEN GOTO 310 

DELWBD = (COSTHETl + COSTHET2 ) / H: GOTO 311 

310 DELWBD =01: 'No downwash contributn. from bound 

' vortex 

311 ' Dummy statement space 

COSl = NUMl / Rl: COS2 = NUM2 / R2 
DELWLH = (11 + COSl) / DENI 
DELWRH = (11 + COS2) / DEN2 

A1(K1, K2) = (DELWLH + DELWBD - DELWRH) / (81 * PI) 
NEXT K2 
NEXT K1 



' ** Similar logic for LH Wing Panel Aero. Influence 
' coeffs. 



315 FOR K1 = 1 TO KMAX: ' K1 Sta. is at the Control Point 
FOR K2 = 1 TO KMAX: ' K2 Sta. is at the Vortex Station 
NUM2 = X3(K1) - X1(K2): NUMl = X3(K1) - X2(K2) 

NUM2 = FC * NUM2: NUMl = FC * NUMl 
' Bring in Prandtl-Glauert factor 
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DEN2 = Y3(K1) + Y1(K2): DENI = Y3(K1) + Y2(K2) 

R1 = SQR(NUM1 " 2 + DENI " 2): R2 = SQR(NUM2"2+DEN2"2 ) 

' Find trig, functns for orthogonal transformation on 
' swept bound vortex 

SINSWP = -SIN(SWPM(K2) ) : COSSWP = COS ( SWPM ( K2 ) ) 

H = NUMl * COSSWP - DENI * SINSWP 

YIROT = NUMl * SINSWP + DENI * COSSWP 

Y2ROT = NUM2 * SINSWP + DEN2 * COSSWP 

COSTHETl = YIROT / R1 ; C0STHET2 = -Y2ROT / R2 

' Logic Check to avoid division by zero 

IF (ABS(H)) <= .001 THEN GOTO 320 

DELWBD = (COSTHETl + COSTHET2) / H: GOTO 321 

320 DELWBD =01: 'No downwash contributn. from bound 

' vortex 

321 ' Dummy statement space 

COSl = NUMl / Rl: COS2 = NUM2 / R2 
DELWLH = (1! + COSl) / DENI 

DELWRH = (11 + COS2) / DEN2 

A2(K1, K2) = (DELWLH + DELWBD - DELWRH) / (8! * PI) 
ASYM(K1, K2) = A1(K1, K2 ) + A2(K1, K2) 

AANT(K1, K2) = A1(K1, K2 ) - A2(K1, K2) 

NEXT K2 
NEXT K1 

400 ' Determine Struct. Infl. Matrix, based on test of Q > 0 
' Note: Q = 0 psi case is rigid wing situation 
IF Q > 01 THEN GOSUB 3000 

' Skip the option for Sym. or Anti-Symm. Soln. 

'GOTO 800: 'Symm. solution branch 

700 ' Select Symmetric or Antisymmetric Solution 
INPUT "Symm. or AntiSymm. Prob. , S/A"; P$ 

710 IF P$ = "S" THEN GOTO 800 
IF P$ = "A" THEN GOTO 900 
GOTO 1000 

800 ' Find Additional Loading 

' GOTO 805: ' Branch to find effect of wing washout 

' Introduce Additional type of alpha 

INPUT "ENTER AOA IN RADIANS, AOA =" , AOA 

FOR K = 1 TO KMAX: ALPHA(K) = AOA: NEXT K 

'GOTO 805 

'PRINT "Elastic wing. Alpha = ", ALPHA(l) 

WRITE #1, ALPHA(l) 

/ 

' Introduce alpha due to built-in geometric twist 
' 0.0 deg. at root, linearly to -2.5 deg. at tip 
WASH = -2.5 * PI / 180! 

FOR I = 1 TO M 

ETA = (I - .5) / M: DELWASH = ETA * WASH 
FOR J = 1 TO N 
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K=(I-1)*N+J: ' Create panel numbering 

ALPHAl(K) = DELWASH: NEXT J 
NEXT I 

'GOTO 805 
/ 

' Introduce alpha (rad.) due to eng. /prop assbly. 

' dead-wgt . 

DATA -0.0012, -0.0042, -0.0084, -0.0138, -0.0189 
DATA -0.0240, -0.0318, -0.0366, -0.0366, -0.0366 
' Read dead-wgt. alphas one at a time. . 

FOR I = 1 TO M 
READ DELALPH 
FOR J = 1 TO N 

K=(I-1) *N+J: ' Create panel numbering 

ALPHA2(K) = NZ * DELALPH: NEXT J 

NEXT I 
/ 

' Introduce alpha (rad. ) due to camber 
DATA -.0525, -.0240, 0.0, .0420, .0730 

' Read in alphas due to camber at the root (WS = 0) 

FOR I = 1 TO N 
READ DELCAM(I) 

NEXT I 

' Calculate the fifty slopes (one per control point) as 
• airfoil varies from 14% to 12% thickness from root to 
' tip 

FOR I = 1 TO M 
FOR J = 1 TO N 
K=(I-1) *N+J 

ALPHA3(K) = DELCAM(J) * (1! - (I - .5) / (7 *M)) 

NEXT J 
NEXT I 

9 

' Introduce alpha (rad. ) due to aileron float angle 
' Calculate aileron float angle (rad) based on velocity 
' (kts) 

FLAIL = -(1.4548 - .0090025 * VEL + .0000444 * VEL " 2) 

* .01745 

FOR I = 1 TO KMAX 
ALPHA4(K) =0!: NEXT I 
ALPHA4(39) = .2 * FLAIL 

ALPHA4(44) = .2 * FLAIL: ALPHA4(49) = .2 * FLAIL 
ALPHA4(35) = .35 * FLAIL 

ALPHA4(40) = FLAIL: ALPHA4(45) = FLAIL: ALPHA4 ( 50 )=FLAIL 

805 ' Combine angles by choice for net alpha distrib. 

INPUT "Enter type of run desired: 0 = Additional only, 

1= Washout only, 2 = Dead wt only, 3 = Camber only, 

4= Aileron float only, 5 = Complete solution, TYPE =", T 
WRITE #1, T 
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FOR I = 1 TO KMAX 
IF (T = 0) THEN 

ALPHA(I) = ALPHA(I); 'Additional loading selection 
ELSEIF (T = 1) THEN 

ALPHA(I) = ALPHAl(I): ' Washout alpha selection 

ELSEIF (T = 2) THEN 

ALPHA(I) = ALPHA2(I): ' Dead-wgt. induced alpha 

' selection 



ELSEIF (T = 3) THEN 

ALPHA(I) = ALPHA3(I): ' Camber alpha section 

ELSEIF (T = 4) THEN 

ALPHA(I) = ALPHA4(I): ' Aileron float selection 

ELSEIF (T = 5) THEN 

ALPHA(I) = ALPHA(I) + ALPHAl ( I ) + ALPHA2 ( I ) + ALPHAS ( I ) 
+ ALPHA4 ( I ) : ' Complete wing selection 

END IF 

XIN(I) = ALPHA(I) 

FOR J = 1 TO KMAX 

A(I, J) = ASYM(I, J) - S(I, J) * Q(I): NEXT J 
NEXT I 



' Find Spanload Soln. L/Q from XIN(KMAX) 

GOSUB 2000 

' S/R returns XIN(I) as L/Q vector of length KMAX 
LIFT =01: MOMENT =01: FOR I = 1 TO KMAX 
LIFT = LIFT + XIN(I) 

MOMENT = MOMENT + XIN(I) * (XMAC25 - .5 * (X1(I) + 

X2 ( I ) ) ) : NEXT I 

' Find Lift & Moment Coefficient 

CL = LIFT * 21 * DELB / S 
/ 

CM = MOMENT * 21 * DELB / (S * MAC) 

NP = .25 - (CM / CL) 

CLTOT = CLREQD - CL 

PRINT "CL ="; CL; ", CM ="; CM; ", CLTOT ="; CLTOT. 

'PRINT "Neut. Pt. at (% Cmac)"; NP: PRINT 
GOSUB 1100: ' Print Header 
CAVE = .5 * CR * (11 + TR) 

FOR I = 1 TO M: ' Sum up L/Q at eta value to get CLC 

ETA = (I - .5) / M: CLC = 01 

FOR J=lTON:K= (I-l) *N+J 

CLC = CLC + XIN(K) 

NEXT J 

' Find struct, twist due to airloads at front panel for 
' sta. "I" 

TWIST =01:J2=(I-1)*N+1 
FOR K = 1 TO KMAX 

TWIST = TWIST + S(J2, K) * Q(I) * XIN(K) : ' Struct twist 

' due to airload 
NEXT K 
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C3 = CR * (11 - (CONSTl * (I-.5)/M)): ' Chord at 

' control point sta. 

CLSECT = CLC / C3 
CLC = CLC / (CL * CAVE) 

GOSUB 1210: ' Print eta, cLc/CLCave, cL(sect) & struct. 

' twist 

NEXT I 
GOTO 990 

900 ' Anti-Symmetric Solution Branch 

' Set up alpha for roll damping due to pB/2V = HELIX 
INPUT "Enter run desired: 0=Roll helix, 
l=Aileron deflection, RUN=", R 
' Calculate available Pb/2V (rad) at specified velocity 
' (kts) 

HELIX = (.21009 - .0008744 * VEL + .000001 * VEL " 2) 

IF (R = 1) GOTO 905 
FOR I = 1 TO M 
ALPDAMP = (I - .5) / M 

' Multiply alpha for roll damping by available roll 
' helix angle (rad) at given velocity (kts) 

ALPDAMP = ALPDAMP * HELIX 
FOR J = 1 TO N 

K= (I-l) *N+J: XIN(K) = ALPDAMP: NEXT J 
NEXT I 
GOTO 908 

905 ' Use following statement for aileron control 
' effectiveness 

FOR I = 1 TO KMAX: XIN(I) =0!: NEXT I 

' Calculate aileron deflection (rad) available at given 
' velocity (kts) 

DELAIL = -(90.30699 - .65754 * VEL + .001771 * VEL " 2 - 
.0000016 * VEL " 3) * .01745 
' Deduct float angle from available control wheel 
' aileron deflection 

FLAIL = -(1.4548 - .0090025 * VEL + .0000444 * VEL " 2) 

* .01745 

DELAIL = DELAIL - FLAIL 

XIN(39) = .2 * DELAIL: XIN(44) = .2 * DELAIL 
XIN(49) = .2 * DELAIL 
XIN(35) = .35 * DELAIL 

XIN(40) = DELAIL: XIN(45) = DELAIL: XIN(50) = DELAIL 

908 ' Solve Spanload Problem 

FOR I = 1 TO KMAX: FOR J = 1 TO KMAX 

A(I, J) = AANT(I, J) - S(I, J) * Q(I): NEXT J 

NEXT I 

' Use S/R ELU/SLVB to solve anti-symm L/Q 
GOSUB 2000 

ROLL =01: FOR I = 1 TO KMAX 
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ROLL = ROLL + XIN(I) * Y3 ( I ) : NEXT I 
CROLL = -ROLL * 2! * DELB / (S * B) 

IF (R = 0) THEN 

PRINT "Roll Damping Deriv. CLP ="; CROLL / HELIX 
ELSE 

PRINT "Aileron Control Deriv., CLa =", CROLL / DELAIL 
END IF 

910 'IF Q > 0 THEN GOTO 915 
'CROLLQO = CROLL 

915 ' Dummy Skip Statement 
'EROLL = CROLL / CROLLQO 
'GOSUB 1200s ' Print results 

' Branch for cLc,CLsect, Twist 
GOSUB 1100s ' Print Header 

CAVE = .5 * CR * (1! + TR) 

FOR I = 1 TO M 'Sum up L/Q at ETA to get CLC 
ETA = (I - .5) / Ms CLC = 0! 

FOR J=lTON: K= (I-l) *N+J 
CLC = CLC + XIN(K) 

NEXT J 

' Find struct, twist due to airloads at front panel for 
' sta. "I" 

TWIST =0lsJ2=(I-l)*N+l 
FOR K = 1 TO KMAX 

TWIST = TWIST + S(J2, K) * Q(I) * XIN(K) 

' Struct twist due to airload 
NEXT K 

C3 = CR * (1! - (CONSTl *(I-.5)/M))s ' Chord at 

' control point sta. 

CLSECT = CLC / C3 

'CLC = CLC / (CL * CAVE) 

GOSUB 1210: ' Print eta, cLc, cL(sect) & struct, twist 
NEXT I 

'GOSUB 1100 

'CAVE = .5 * CR * (1! + TR) 

'FOR I = 1 TO M: ' Defaulted for N = 1 (Chrdwse sta. ) 

'ETA = (I - .5) / M 

'CLC = XIN(I) / CAVE: GOSUB 1200 

'NEXT I 



990 ' Q = Q + i: 
'NEXT MM 
'GOTO 700 



1000 END 
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1100 ' S/R to Print header, whole selection shown 

'PRINT TAB(3); "Mach No."; TAB (14); "CLalpha"; TAB (25) 
' "CmCl"; TAB(35); "N.P."; TAB(45); "Clp" 

'PRINT TAB(5); "Q, psi"; TAB (14); "CLalpha"; TAB (25); 

' "CmCl"; TAB(35); " N.P."; TAB(45); 

' "Elift"; TAB (55); "Del-NP" 

'PRINT TAB(IO); "Roll Control Evaluation, Mid Ail. at 
' 0.05" 

'PRINT TAB(5); "Q, psi"; TAB(15); "Croll"; TAB(25); 

' "Eroll" 

PRINT 

PRINT TAB(6); "2Y/B"; TAB (13); "Clc/CLcave" ; TAB (26); 

' "Cl"; TAB(35); "Twist"; TAB(45); "Degrees" 

RETURN 

1200 ' S/R to Print results, whole selection shown 

'PRINT USING "#####.####"; Q; CL; CM; NP; ELIFT; DELNP 
'PRINT USING "#####.####"; Q; CROLL; EROLL 
1210 PRINT USING "#####.####"; ETA; CLC; CLSECT; TWIST; 

TWIST / .01745 

WRITE #1, ETA, CLC, CLSECT, TWIST, TWIST / .01745 
RETURN 

2000 ' S/R ELU 

' Tri-Diagonalizes the input matrix A(KMAX,KMAX) 

' Input Column vector, XIN(KMAX), gets replaced by 
' the output which is returned to the main program 
NMl = KMAX - 1 

FOR K = 1 TO NMl: KPl = K + 1 
FOR I = KPl TO KMAX 

G = -A(I, K) / A(K, K): A(I, K) = G 
FOR J = KPl TO KMAX 

A(I, J) = A(I, J) + G * A(K, J): NEXT J 
NEXT I 
NEXT K 

GOSUB 2050: ' Use S/R SLVB for next step 

RETURN 

2050 ' S/R SLVB 

' Solves the Tri-Diagonalized matrix [A] obtained 

' from ELU by back substitution 

NMl = KMAX - 1: NPl = KMAX + 1 

FOR K = 1 TO NMl 

KPl = K + 1 

FOR I = KPl TO KMAX 

XIN(I) = XIN(I) + A(I, K) * XIN(K) 

NEXT I 
NEXT K 

XIN(KMAX) = XIN(KMAX) / A(KMAX, KMAX) 

FOR K = 2 TO KMAX 
I = NPl - K 
J1 = I + 1 
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FOR J = J1 TO POIAX 

XIN(I) = XIN(I) - A(I, J) * XIN(J) 

NEXT J 

XIN(I) = XIN(I) / A(I, I) 

NEXT K 
RETURN 

3000 'S/R to Determine [S] Structural Infl. Coeff, matrix 
/ 

3010 ' Generate [My] and [Mx] matrices, order M x KMAX 
' First get diagonal elements 
K = 1: FOR I = 1 TO M: FOR J = 1 TO N 
MY(I, K) = .5 * (XEA(I) - X1(K) - .75 * DELB * 
TAN(SWP(K) ) ) 

MX(I, K) = .25 

K = K + 1: NEXT J 

NEXT I 

' Now get Upper triangular elements 
VALUE = 1 ! : FOR L = 1 TO M - 1 
K=L*N+1; K2=M-L 
FOR I = 1 TO K2; FOR J = 1 TO N 

MY(I, K) = XEA(I) - X1(K) - .5 * DELB * TAN(SWP(K)) 

MX(I, K) = VALUE: K = K + 1: NEXT J 

NEXT I 

VALUE = VALUE + 1 ! : NEXT L 
' Test circuit 

' PRINT "Print out of MX(I,J)" 

' FOR I = 1 TO M: FOR J = 1 TO KMAX 
' PRINT USING "##.##"; MX(I, J); 

' NEXT J 

' PRINT : NEXT I 

3020 ' Generate [DELA(I,J)], order KMAX x M 
K = 1 

FOR I = 1 TO- M: FOR J = 1 TO N 
DELA(K, I) = 1! 

K = K + 1.: NEXT J 
NEXT I 

' Test circuit 

' PRINT "Print out of DELA(I,J)" 

' FOR I = 1 TO ICMAX; FOR J = 1 TO M 
' PRINT USING "##.##"; DELA(I, J); 

' NEXT J 

' PRINT : NEXT I 



3030 ' Generate [U(I,J)] square matrix, order M x M 
FOR I = 1 TO M: FOR J = 1 TO M 

IF I = J THEN U(I, J) = .5 

IF I < J THEN U(I, J) = 0! 

IF I > J THEN U(I, J) = 1! 
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NEXT J 
NEXT I 

' Test circuit 

' FOR I = 1 TO M: FOR J = 1 TO M 
' PRINT USING "##.##"; U(I, J); 

' NEXT J 

' PRINT : NEXT I 

3040 ' Generate UGJ(I,J) & UEI(I,J) matrices, order M x M 
' Embed DELB^2 type of constants 
DELB2 = DELB " 2: TANDELB2 = TANEA * DELB2 
FOR I = 1 TO M: FOR J = 1 TO M 

' DELB2 =1!: GJ(J) =6! -J : ' for test circuit 
UGJ(I, J) = DELB2 * U(I, J) / GJ(J) 

UEI(I, J) = TANDELB2 * U(I, J) / EI(J): NEXT J 
NEXT I = TANDELB2 * U(I, J) / EI(J): NEXT J 

'Test Circuit 

' PRINT "Sample of UGJ(I,J)" 

' FOR I = 1 TO M: FOR J = 1 TO M 
' PRINT USING ••##.##"; UGJ(I, J); 

' NEXT J 
'PRINT : NEXT I 

3050 ' Start gathering terms for [S] matrix assembly 
DELSIN = DELB * SINEA: DELCOS = DELB * COSEA 



FOR I = 


1 TO M: FOR J = 1 


TO 


KMAX 




F(I, J) 


= 0! : G(I, J) =01 








FOR K = 


1 TO M 








F(I, J) 


= F(I, J) + UGJ(I, 
DELSIN * MX(K, J) ) 


K) 


* (COSEA * 


MY(K, J) + 


G(I, J) 


= G(I, J) + UEI(I, 


K) 


* (SINEA * 


MY(K, J) - 



DELCOS * MX(K, J) ) 

NEXT K 
NEXT J 
NEXT I 

3060 ' Form the [S] matrix, order KMAX x KMAX 
FOR I = 1 TO KMAX: FOR J = 1 TO KMAX 
S(I, J) = 01 
FOR K = 1 TO M 

S(I, J) = S(I, J) + DELA(I, K) * (F(K, J) + G(K, J) ) : 
NEXT K 
NEXT J 
NEXT I 

' [S] matrix is now available 
RETURN 
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APPENDIX B 



P-3 LEADING EDGE FAILURE ANALYSIS 
TWO DIMENSIONAL PANEL METHOD 

PROGRAM PANEL 

.THIS FORTRAN PROGRAM ITERATIVELY COMPUTES THE PRESSURE 
AND FORCE DISTRIBUTION AROUND THE P-3 AIRFOIL AT 24 
SELECTED WING STATIONS AND WRITES A LOAD FILE FOR USE 
IN MSC PAL2 FINITE ELEMENT ANALYSIS SOFTWARE. INPUT 
CONSISTS OF THE AIRFOIL COORDINATES (X,Y) FROM A FILE 
CALLED P3.DAT, THE 24 WING STATION LOCATIONS IN 
INCHES AND THE NODE NUMBERS FROM THE FINITE ELEMENT 
MODEL AT WHICH FORCES AND CONSTRAINTS WILL BE APPLIED. 
OUTPUT CONSISTS OF THE PAL2 LOAD FILE (LE.LD), A 
TABULAR FILE IDENTIFYING THE COEFFICIENT OF PRESSURE, 
DYNAMIC PRESSURE AND FORCES AT A GIVEN PANEL MID POINT 
ON THE LEADING EDGE OF THE AIRFOIL. SECTION LIFT 
COEFFICIENT (Cl) IS MATCHED TO AN INPUT Cl REQUIREMENT 
(GENERATED BY AN AEROELASTIC SPAN LOAD ANALYSIS 
PROGRAM) THROUGH ITERATION OF ANGLE OF ATTACK. 

PARAMETER(N=53,M=N+1,PI=3. 14159265,RHO=. 0023769, 

+ PAMB=14.7) 

REAL X(M) ,Y(M) ,THETA(M) ,BETA(M,M) ,R(M,M) ,XM(M) , YM(M) , 

+ A(100, 101) ,AN(M,M) ,BN(M,M) , AT(M,M) ,BT(M,M) , 

+ SUMBN(M,M) ,SUMB1,SUMB2,SUMA(M) ,SUMB(M) ,VTAN(M) , 

+ Q(M) ,CP(M) ,P(M) ,NUM(M) ,DEN(M) , ALPHA, V,VEL,AOA, 

+ CY(M) ,CX(M) ,CM(M) ,FY(M) ,FX(M) , VTOT (N) , CYT, CXT, CMT 

+ FYT,FXT,CD,CL,CLRQD(100) ,WS(100) ,XML(M) ,YML(M) , 

+ LOADM(M) , INTWIST, OUTWIST, TWIST, XTOPDIS,XBOTDIS, 

+ XDISTOP(48) ,XDISHIN(24) , ZDISTOP( 48 ) , ZDISHIN( 24 ) 

INTEGER K,I, J, NODE ( 100, 100) ,NTOP( 100) ,NHIN( 100) 

OPEN (UNIT=51,FILE='P3. DAT ' ,STATUS=' UNKNOWN ' ) 

OPEN (UNIT=5 7, FILE=' DUMP. DAT ' ,STATUS=' UNKNOWN ' ) 

OPEN (UNIT=58,FILE='WS. DAT ' ,STATUS=' UNKNOWN' ) 

OPEN ( UNIT=6 0, FILE= ' LE. LD ' ,STATUS=' UNKNOWN ' ) 

OPEN (UNIT=61,FILE=' TOP. DAT ' ,STATUS=' UNKNOWN' ) 
OPEN(UNIT=62,FILE='HINGE.DAT' ,STATUS='UNKNOWN' ) 

OPEN (UNIT=6 3, FILE=' NODESa.DAT' , STATUS = 'UNKNOWN ' ) 

OPEN (UNIT=64 , FILE= ' NODESb . DAT ' , STATUS= 'UNKNOWN ' ) 

OPEN (UNIT=6 5, FILE=' NODESc.DAT' , STATUS= ' UNKNOWN ' ) 
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READ IN AIRFOIL DATA 



DO 10 K=1,N 

READ(51,500) X(K),Y(K) 

10 CONTINUE 
500 FORMAT (2F10. 6) 

X(N+1)=X(1) 

Y(N+1)=Y(1) 

C READ IN FINITE ELEMENT NODE NUMBERS FOR LATER USE 

C (READ AS NODE (PANEL NO., RIB)) 

DO 12 1=17,36 

READ (63, 575) (NODE ( I , J ) , J=1 , 8 ) 

READ(64,575) ( NODE ( I , J ) , J=9 , 16 ) 

READ(65,575) (NODE ( I , J ) , J=17 , 24 ) 

12 CONTINUE 
575 FORMAT(9I5) 

C WRITE HEADER FOR FINITE ELEMENT LOAD FILE 

WRITE (60, 590) 

C REQUEST AND RECEIVE AIRSPEED INPUT 

PRINT*, 'ENTER VELOCITY IN KNOTS (XXX. XX) ' 

READ*,VEL 

V = VEL*1.6878 

C REQUEST AND RECEIVE FRONT SPAR TWIST ANGLES 

PRINT*, 'ENTER TWIST AT 2Y/B=.35 ( rad) ( .XXXX) ' 
READ*,TWIST35 

PRINT*, 'ENTER TWIST AT 2Y/B=.45 ( rad) (. XXXX) ' 
READ*,TWIST45 

PRINT*, 'ENTER TWIST AT 2Y/B=.55 ( rad) ( .XXXX) ' 
READ*,TWIST55 



C CALCULATE CHORD LENGTH OF P-3 WING AT SELECTED WING 

C STATIONS 



C SET UP OUTER LOOP TO ITERATIVELY CALCULATE PRESSURES 

C FOR A SERIES OF WING STATIONS 
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INPUT WING STATIONS, Cl AND LOAD MULTIPLIER FROM FILE 

WS.DAT AND CALCULATE CHORD LENGTH AT THE WING STATION 

502 FORMAT(F7.1,F8.4,F6.2) 

DO 15 L=l,24 

READ(58,502) WS (L) ,CLRQD(L) ,LOADM(L) 

PRINT* , ' ITERATION= ' , L 
PRINT*, 'WS=' ,WS(L) 

PRINT* , ' CLRQD= ' , CLRQD ( L ) 

CHORD = 227.*(1-(0.6*(2.*WS(L)/1188. ) ) ) 

PRINT* , ' CHORD= ' , CHORD 

CALCULATE THE THICKNESS FRACTION AT A GIVEN WING 

STATION (USE TO LINEARLY VARY AIRFOIL THICKNESS WITH 
WING STATION) 

DTHICK = l-( (1-(12./14. ) )*(2.*WS(L)/1188. ) ) 

PRINT*, 'THICKNESS RATIO =', DTHICK 

SCALE THE THICKNES OF THE AIRFOIL AT EACH WING STATION 

BY THE THICKNESS FRACTION 

DO 16 K=1,N 

Y(K)=Y(K)*DTHICK 
16 CONTINUE 

C SET STARTING ALPHA 

ALPHA = (PI/180. ) 

C SET COUNTER START FOR AOA ITERATION 

KOUNT = 1 

C BEGIN ITERATION OF ALPHA TO MATCH CLREQD = CL 

1100 CONTINUE 

C FIND NORMALIZED MID-POINT OF PANEL 

DO 30 1=1, N-1 

XM(I)=0.5*(X(I)+X(I+1) ) 

YM(I)=0.5*(Y(I)+Y(I+1) ) 

C FIND MID-POINT OF PANEL IN INCHES 

XML ( I ) =XM ( I ) *CHORD 
YML(I)=YM(I)*CHORD 
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CALCULATE THETA (ANGLE BETWEEN PANEL AND CHORD LINE) 



NUM(I)=Y(I+1)-Y(I) 

DEN(I)=X(I+1)-X(I) 

THETA(I)=ATAN2(NUM(I) ,DEN(I) ) 

C CALCULATE R (DISTANCE BETWEEN PANEL MID-POINTS) 

DO 20 J=1,N 

R(I, J)=SQRT( (XM(I)-X( J) )**2+(YM(I)-Y( J) )**2) 

C CALCULATE BETA (ANGLE AT MID POINT OF PANEL 1, 

C ENCLOSING END POINTS OF PANEL J) 

BETA(I, J)=ATAN2( ( ( YM( I )-Y( J+1 ) ) * (XM( I ) -X( J ) ) 

# -( (XM(I)-X(J+1) ))*(YM(I)-Y(J))), 

# ( (XM(I)-X( J+1) )*(XM(I)-X( J) ) 

# +(YM(I)-Y( J+1) )*(YM(I)-Y( J) ) ) ) 



20 CONTINUE 
30 CONTINUE 

C CALCULATE AUGMENTED "A" MATRIX 



SUMB1=0.0 
SUMB2=0.0 
DO 50 1=1, N-1 
SUMBN ( I , N ) =0 . 0 
DO 40 J=1,N-1 

IF (I.EQ.J) THEN 
AN(I, J)=0.5 
AT(I,J)=0.0 
BN(I, J)=0.0 
BT(I, J)=0.5 
ELSE 

AN(I, J)=(SIN(THETA(I)-THETA( J) ) *LOG( R( I , J+1 ) /R( I , J) ) 

# +COS(THETA(I)-THETA( J) )*BETA(I, J) )/(2*PI) 

BN(I, J)=(COS(THETA(I)-THETA( J) ) *LOG( R( I , J+1 ) /R( I , J) ) 

# -SIN(THETA(I)-THETA( J) )*BETA(I, J) )/(2*PI) 
AT(I, J)=-BN(I, J) 

BT(I, J)=AN(I, J) 

END IF 

SUMBN(I,N)=SUMBN(I,N) + BN(I,J) 

IF (I.EQ.l) THEN 

SUMB1=SUMB1 + BT(1,J) 

ELSEIF (I.EQ.N-1) THEN 
SUMB2=SUMB2 + BT(N-1,J) 

END IF 
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A(I,J) = AN(I,J) 

A(I,N) = SUMBN(I,N) 

A(N,J) = AT(1,J) + AT(N-1,J) 

A(N,N) = SUMBl + SUMB2 
40 CONTINUE 
50 CONTINUE 

DO 60 1=1, N-1 

A(I,N+1) = -V*SIN(ALPHA-THETA(I) ) 

60 CONTINUE 

A ( N , N+ 1 ) =-V* ( COS ( ALPHA-THETA ( 1 ) ) + 

+ COS ( ALPHA-THETA ( N- 1 ) ) ) 

CALL GAUSS (N, A) 

GAMMA=0 . 0 
DO 80 1=1, N 
Q(I)=A(I,N+1) 

GAMMA=A(N,N+1) 

80 CONTINUE 

DO 100 1=1, N-1 
SUMA(I)=0.0 
SUMB(I)=0.0 
DO 90 J=1,N-1 

SUMA( I )=SUMA( I ) +Q ( J ) *AT ( I , J ) 

SUMB ( I ) =SUMB ( I ) +BT ( I , J ) 

VTOT ( I ) =SUMA ( I ) +GAMMA*SUMB ( I ) +V*COS ( ALPHA-THETA ( I ) ) 
VTAN(I)=VTOT(I)/V 
CP(I)=1.0-(VTAN(I) )**2.0 
P(I)=C(CP(I)*0.5*RHO*(V**2. ) )/144. ) 

90 CONTINUE 
100 CONTINUE 

C CALCULATE FORCE COEFFICIENTS (CX,CY) AND FORCES (FX= 

C CHORDWISE, FY=NORMAL TO CHORD) AT PANEL MID POINT 



C 



CYT=0 

CXT=0 

CMT=0 

FXT=0 

FYT=0 

DO 110 1=1, N-1 

CY(I)=-1.0*CP(I)*(X(I+1)-X(I) ) 
CX(I)=CP(I)*(Y(I+1)-Y(I) ) 

MULTIPLY BY CHORD FOR FULL SCALE SOLUTION 



FY(I)=CHORD*CY( I)*0.5*RHO*(V**2. )/144 
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FX(I)=CHORD*CX(I)*0.5*RHO*(V**2. )/144 
CM(I)=CP(I)*( (X(I+1)-X(I) )*XM(I)+ 

+ (Y(I+1)-Y(I) )*YM(I) ) 

CYT=CYT+CY(I) 

CXT=CXT+CX ( I ) 

FYT=FYT+FY(I) 

FXT=FXT+FX(I) 

110 CONTINUE 

C GET FORCE TOTALS FORWARD OF FRONT SPAR (.15C) 

FYNT=0 . 0 
FXNT=0 . 0 
DO 115 1=18,35 
FYNT=FYNT+FY ( I ) 

FXNT=FXNT+FX(I) 

115 CONTINUE 

CALCULATE APPROXIMATE TOTAL FORCES ON COMPLETE LE 
SECTION (BETWEEN NACELLES) BY TAKING FORCES AT 
MID POINT * 92 INCHES 

IF (L.EQ.7) THEN 
FYTOT=FYNT*92 
FXTOT=FXNT*92 
END IF 

C PERFORM AXIS ROTATION TO FIND CL 

CL=CYT* COS ( ALPHA ) -CXT* S IN ( ALPHA ) 

C PERFORM ITERATION OF ANGLE OF ATTACK 

DIFF=.0001 
CLDIFF=CLRQD ( L ) -CL 
KILL=50 
DELTA=0.1745 
IF (KOUNT.GT.KILL) THEN 
PRINT*, 'COUNTER EXCEEDED' 

GOTO 111 

ELSEIF (CLDIFF.GT.DIFF) THEN 

PRINT' (I3,2F10.6) ' , KOUNT , ALPHA , CL 
KOUNT=KOUNT+l 

ALPHA=ALPHA+ ( DELTA* ABS ( CLDIFF ) ) 

GOTO 1100 

ELSEIF (CLDIFF.lt. 0.0) THEN 

PRINT' (I3,2F10.6) ', KOUNT , ALPHA , CL 
KOUNT=KOUNT+l 

ALPHA=ALPHA-(DELTA*ABS( CLDIFF) ) 

GOTO 1100 
ELSE 
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' (4F10.5) ' , ALPHA, CL, FXNT^FYNT 
(57,570) VEL,WS(L) , ALPHA, CLRQD 
(57,560) 



PRINT 

WRITE \ o ! , o !\j } 
WRITE (57, 560) 
ENDIF 



(L) ,CL,FXNT,FYNT 



C WRITE PRESSURES AND FORCES TO OUTPUT FILES 



DO 120 1=17,36 
WRITE (57, 580) 

+ 

120 CONTINUE 



I,XM(I),XML(I),YML(I),CP(I),P(I), 
FX(I)*LOADM(L) ,FY(I)*LOADM(L) 



C WRITE LOADING FILE 

APPLICATION ( PAL2 ) . 



(LE.LD) FOR FINITE ELEMENT 



FX( 


17 


)=FX( 


17 


)*0. 


5 




FY( 


17 


)=FY( 


17 


)*0. 


5 




FX( 


23 


)=FX( 


23 


) +0 . 


5*FX 


(24) 


FY( 


23 


)=FY( 


23 


)+0. 


5*FY 


(24) 


FX( 


25 


)=FX( 


25 


)+0. 


5*FX 


(24) 


FY( 


25 


)=FY( 


25 


)+0. 


5*FY 


(24) 


FX( 


26 


).=FX( 


26 


)+FX 


:(27) 




FY( 


26 


)=FY( 


26 


)+FY 


(27) 




FX( 


28 


)=FX( 


28 


)+0. 


5*FX 


(29) 


FY( 


28 


)=FY( 


28 


)+0. 


5*FY 


(29) 


FX( 


30 


)=FX( 


30 


) ■*■0 c 


5*FX 


(29) 


FY( 


30 


)=FY( 


30 


) +0 . 


5*FY 


(29) 


FX( 


36 


)=FX( 


36 


)*0. 


5 




FY( 


36 


)=FY( 


36 


)*0. 


5 




DO 


125 1=17, 


36 







J=I 

IF (I.LE.23) THEN 

WRITE(60,591) NODE(J,L) 

NODE ( J , L ) 
EQ.25) THEN 

NODE ( J , L ) 
NODE( J,L) 
26) THEN 

NODE ( J , L ) 
NODE ( J , L ) 
(I.EQ.28) THEN 

NODE ( J , L ) 
NODE ( J , L ) 
(I.GE.30) THEN 

NODE ( J , L ) 
NODE( J,L) 



ELSEIF (I 

WRITE(60,591) 

h 

ELSEIF (I.EQ 
WRITE (60, 591) 

h 

ELSEIF 

WRITE (60, 591) 

h 

ELSEIF 

WRITE(60,591) 



,FX(I)*LOADM(L) , 
,FY(I)*LOADM(L) 

,FX(I)*LOADM(L) , 
,FY(I)*LOADM(L) 

,FX(I)*LOADM(L) , 
,FY(I)*LOADM(L) 

,FX(I)*LOADM(L) , 
,FY(I)*LOADM(L) 

,FX(I)*LOADM(L) , 
,FY(I)*LOADM(L) 
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ENDIF 

CONTINUE 
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560 FORMAT(2X, 'PANEL' , 3X, 'XM( I ) ' , 4X, 'XML(I) ' ,4X, 'YML(I) ' , 
+ 5X, 'CP(I) ',5X, 'P(I) ',6X, 'FX(I) ',5X, 'FY(I) ') 

570 FORMAT(//,2X, 'VEL= ' , FIO . 1 , / , 2X, 'WS= ' , FIO . 3 , / , 2X, 

+ 'ALPHA=',F10.5,/,2X, ' CLRQD= ' , FIO . 5 , / , 2X, 'CL=' , 

+ F10.5,/,2X, 'FXNT=' ,F10.5,/,2X, ' FYNT= ' , FIO . 5 , / ) 

580 FORMAT(2X,I3,7F10.5) 

590 FORMAT( IX, 'FORCES AND MOMENTS APPLIED O') 

591 FORMAT (IX, ' FX ' , 17 , FIO . 5 , / , IX, ' FZ ' , 17 , FIO . 5 ) 

111 CONTINUE 

15 CONTINUE 

C CALCULATE TWIST DISPLACEMENTS OF FINITE ELEMENT NODES 

INTWIST=0 . 82155* ( TWIST45-TWIST35 ) +TWIST35 

OUTWIST=0.89371*(TWIST55-TWIST45)+TWIST45 

TWIST=OUTWIST-INTWIST 

XTOPDIS=TWIST*8 . 8 

XBOTDIS= (-1.0) *XTOPDIS 

DO 130 1=1,24 
J=I-1 

XDISTOP(2*I-l)=(XTOPDIS/23. )*J 
XDISTOP(2*I)=(XTOPDIS/23. ) * J 
XDISHIN(I)=(XBOTDIS/23)*J 
TT=(TWIST/23. )*J 
XT=(XTOPDIS/23. ) * J 
IF (TT.LE.0625) THEN 
ZDISTOP(2*I)=0.5*TT*XT 
ELSE 

ZDISTOP(2*I)=(-0.5)*(TT-.125)*XT 
END IF 

ZDISTOP( 2*1-1 ) = ( -0 . 5 ) *TT* (XT+1 . 0 ) 

ZDISHIN( I )=0 . 5*TT*XT 
130 CONTINUE 

WRITE (5 7, 55 5) FXTOT,FYTOT 
555 FORMAT(/,2X, 'APPROX TOTAL FX( 92" ) = ' / FIO . 1 , 

+ /,2X, 'APPROX TOTAL FY( 92" ) = ' , FIO . 1 ) 

C READ AND WRITE NODES AND CORRESPONDING BOUNDARY 

C CONDITIONS IN FINITE ELEMENT LOAD 

WRITE (60, 597) 

READ(61,592) (NTOP( I ) , 1=1 , 48 ) 

READ(62,601) (NHIN( I ) , 1=1 , 24 ) 

WRITE (60, 593) (NTOP( I ) , 1=1 , 48 ) 

WRITE (60, 5931) (NTOP(I) ,XDISTOP(I) ,1=1,48) 
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n n non 



ENFORCE Y-AXIS DISPLACEMENT CONSTRAINT ALONG FRONT 

EDGE OF SPAR CAP FLANGE ONLY. (RELEASE BACK EDGE TO 
SIMULATE SINGLE LINE OF SCREW FASTENERS.) 

WRITE (60, 5932) ( NTOP ( 2*1 ) , 1=1 , 24 ) 

WRITE (60, 5933) ( NTOP ( I ) , ZDISTOP ( I ) , 1=1 , 48 ) 

WRITE (60, 594) ( NHIN ( I ) , XDISHIN ( I ) ,1=1,24) 

WRITE(60,599) (NHIN( I ) , 1=1 , 24 ) 

WRITE(60,600) (NHIN( I ) , ZDISHIN( I ) , 1=1 , 24 ) 
WRITE(60,595) (NHIN( I ) , 1=1 , 24 ) 

WRITE(60,596) (NHIN( I ) , 1=1, 24) 

WRITE (60, 598 ) 

592 FORMAT (15) 

601 FORMAT (15) 

593 FORMAT(lX, 'RA' ,I7,1X, '0' ) 

5931 FORMAT(lX, 'TX',I7,1X,F8.4) 

5932 F0RMAT(1X,'TY',I7,1X, '0') 

5933 FORMAT(lX, 'TZ' , I7,1X,F8.4) 

594 FORMAT(lX, 'TX',I7,1X,F8.4) 

599 F0RMAT(1X, 'TY' ,I7,1X, '0^ ) 

600 F0RMAT(1X, 'TZ' ,I7,1X,F8.4) 

595 FORMAT(lX, 'RX',I7,1X, ^0' ) 

596 FORMAT(lX, 'RZ' ,17, IX, '0' ) 

597 FORMAT(lX, ' — BLANK LINE — 

+ IX, 'DISPLACEMENTS APPLIED O') 

598 FORMAT(lX,' — BLANK LINE — ',/, 

+ IX, 'SOLVE' ,/, IX, 'QUIT' ) 



STOP 

END 



SUBROUTINE GAUSS (N, A) 

...THIS SUBROUTINE PERFORMS GAUSS ELIMINATION WITH 
PIVOTING. 

INTEGER PV ’ ! PIVOT INDEX. 

DIMENSION A( 100, 101) 

EPS=1.0 

10 IF (1.0+EPS.GT. 1.0) THEN 
EPS=EPS/2.0 
GOTO 10 
END IF 
EPS=EPS*2.0 
EPS2=EPS*2.0 
DO 1010 1=1, N-1 
PV=I 

DO J=I+1,N 

IF (ABS(A(PV,I) ) .LT.ABS(A(J,I) ) ) PV=J 
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END DO 

IF (PV.EQ.I) GOTO 1050 
DO JC=1,N+1 
TM=A(I, JC) 

A(I, JC)=A(PV, JC) 

A(PV, JC)=TM 
END DO 

1050 IF (A(I,I) .EQ.O) GOTO 1200 ! KICK OUT IF SINGULAR 

DO JR=I+1,N ! ELIMINATION OF BELOW DIAGONAL. 

IF (A(JR,I) .NE.0.0) THEN 
R=A(JR, I)/A(I, I) 

DO KC=I+1,N+1 
TEMP=A( JR,KC) 

A(JR,KC)=A( JR,KC) - R*A(I,KC) 

IF (ABS(A(JR,KC) ) .LT.EPS2*TEMP) A ( JR, KC )=0 . 0 

C IF THE RESULT OF SUBTRACTION IS SMALLER THAN 2 TIMES 

C EPS TIMES THE ORIGINAL VALUE, IT IS SET TO ZERO. 

END DO 
END IF 
1060 END DO 
1010 CONTINUE 

IF (A(N,N) .EQ.O) GOTO 1200 
A(N,N+1)=A(N,N+1)/A(N,N) 

DO NV=N-1,1,-1 ! BEGIN BACKWARD SUBSTITUTION. 

VA=A(NV,N+1) 

DO K=NV+1,N 

VA=VA-A ( NV, K ) * A ( K , N+1 ) 

END DO 

A ( NV, N+1 ) =VA/A ( NV, NV ) 

END DO 
RETURN 

1200 PRINT*, 'MATRIX IS SINGULAR' 

STOP 

END 
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